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ABSTRACT 

Cosmic dust is present in many astrophysical objects, and recent observations across the electromag- 
netic spectrum have revealed that the dust distribution is often strongly three-dimensional. Dust grains 
are effective in absorbing and scattering UV/optical radiation, and re-emit the absorbed energy at infrared 
wavelengths. Understanding the intrinsic properties of these objects, including the dust itself, therefore 
requires 3D dust radiative transfer calculations. Unfortunately, the 3D dust radiative transfer problem is 
non-local and non-linear, which makes it one of the hardest challenges in computational astrophysics. 
Nevertheless, significant progress has been made in the last decade, with an increasing number of codes 
capable of dealing with the complete 3D dust radiative transfer problem. We discuss the complexity of 
this problem, describe the two most successful solution techniques (Ray-Tracing and Monte Carlo), and 
discuss the state of the art in modeling observational data using 3D dust radiative transfer codes. We end 
with an outlook on the bright future of this field. 

Subject headings: scattering - Monte Carlo - ray-tracing - computational astrophysics - numerical algorithms 



1. INTRODUCTION 

Given the dominant role of radiation in astro- 
physics, its transport through a medium is one of the 
most fundamental processes to be considered. Ana- 
lyzing the radiation received from an object not only 
provides us with information about its radiation source 
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but also the medium in between and surrounding the 
object and the observer. 

Among the numerous ways of producing and pro- 
cessing radiation, dust grains mixed in the cosmic gas 
play a special role. They are efficient at absorbing 
and scattering ultraviolet (UV) through near-infrared 
(NIR) photons and then re-radiating the absorbed en- 
ergy in the infrared and sub-millimeter (submm) wave- 
length range. Cosmic dust can be found in many as- 
trophysical objects like the solar system (Hoppe et al. 
2010), comets and meteoroids (Kiippers et al. 2005), 
sub-stellar atmospheres (Harvey et al. 2012), young 
stellar objects (Keller et al. 2008), proto-stellar to 
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proto-planetary disks (Watson et al. 2009), evolved 
stars (Groenewegen et al. 2011), reflection nebulae 
(Castellanos et al. 2011), supernova remnants (Rho 
et al. 2009), molecular clouds (Martel, Urban & Evans 
2012), the interstellar medium (Zhukovska, Gail & 
Trieloff 2008), galaxies (Dunne et al. 2011), active 
galactic nuclei (Haas et al. 2000), and the high-redshift 
universe (Dwek, Galliano & Jones 2007). As dust 
grains modify the radiation field in these objects, they 
should be taken into account for an unbiased analysis 
of their intrinsic properties. Such an analysis requires 
radiative transfer (RT) calculations to be performed. 

Aside from its importance as a tracer, the physical 
and chemical processes related to dust itself are of in- 
terest. They cover its formation, its cycle in galaxies, 
the variation in opacity with chemical composition, its 
growth and destruction processes in cloud cores and 
circumstellar disks to act as building blocks for plan- 
ets, its interaction with magnetic fields and the chem- 
istry on the grain surface. For example, the dust RT is 
important for understanding chemistry in the ISM as 
photo-dissociation rates are strongly dependent on the 
UV radiation field that includes a significant amount of 
photons scattered from dust grains. In this review, we 
will just address physical properties of the dust where 
they enter the RT or the modeling of objects, and refer 
the reader to one of the many published works on dust 
for the aspects mentioned above (e.g., Draine 2003a; 
Henning, Griin & Steinacker 2009; Henning 2010). 

Many dusty objects have been observed at in- 
creasingly higher spatial resolution in the last ten 
years. These observations cover UV/optical/NIR (e.g., 
Hubble, GALEX, ground-based telescopes) and in- 
frared/submm wavelengths (e.g., Spitzer, Akari, Her- 
schel, Planck, WISE, ALMA, JCMT, APEX, IRAM). 
Space telescopes exploring atmospherically absorbed 
wavelength windows, high-resolution interferomet- 
ric data, polarization data, or all-sky maps are just 
a few examples of the rich data set that awaits the 
RT modeler. The information determined by compar- 
ison of dust RT calculations with global and pixel- 
by-pixel resolved spectral energy distributions (SEDs) 
of dusty objects include the properties of the illumi- 
nating sources (stars, accretion disks, integrated star 
formation rate, etc.), the distribution of the dust (disk 
structures, cloud geometries, underlying multi-phase 
nature, etc.), and properties of the dust grains (size, 
shape, and composition). A feature commonly seen 
in high spatial resolution images at all wavelengths 
is the complex nature of the dust density distribution. 



Examples of global 3D geometries include complex 
arm structures in spiral galaxies (Patrikeev et al. 2006; 
Fritz et al. 2012), large scale filaments in star forming 
regions (Andre et al. 2010; Arzoumanian et al. 2011), 
and bow-shocked shells around evolved stars (Cox 
et al. 2012). A prominent example of locally complex 
3D geometries is the known fractal nature of the in- 
terstellar medium (Beech 1987; Falgarone, Phillips & 
Walker 1991). Images illustrating complex local and 
global 3D dust structures are shown in Fig. 1. In ad- 
dition, the illuminating sources of dust have been long 
know to have complex distributions from the com- 
bination of the anisotropic interstellar radiation field 
and local neighboring stars to the stellar distribution 
in galaxies. Both of these issues show that a complete 
3D treatment of RT is inevitable and critical to make 
progress in many fields. 

Among the many computational problems in astro- 
physics, 3D line and dust RT has long been a major 
challenge and often approximated or neglected. While, 
for example, 3D (magneto-)hydrodynamics ((M)HD) 
codes have existed for many years, radiation transport 
is considered to be one of the four Grand Challenges 
in Computational Astrophysics ~. Dust RT is different 
than line RT in that the dust opacities generally do not 
depend on the RT solution itself. 

The reasons to neglect or approximate 3D dust RT 
are manifold. A good portion of the difliculties arises 
from the fact that the underlying physical processes 
combine, in the stationary case, to a non-local and 
non-linear 6D problem. Since the radiation field needs 
to be determined in all directions, at any spatial loca- 
tion, and for each wavelength, the solution vector itself 
comprises three dimensions more than the variables 
in (M)HD problems. The RT problem is non-local 
in space (propagation of the photons within the en- 
tire domain), direction (scattering and absorption/re- 
emission), and wavelength (absorption/re-emission). 
This non-locality makes it difficult to simplify the 
problem by neglecting processes or wavelengths. 
For example, absorption and scattering have roughly 
the same efficiency from UV to NIR, with strongly 
anisotropic scattering (Gordon 2004). Modeling far- 
infrared (FIR) images, a consistent treatment of the 
dust emission requires the RT to be calculated where 
the dust absorption happens, at shorter wavelengths. 



^see, e.g., the Grand Challenge conference series at the Institute for 
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Fig. 1. — The complex and filamentary structures of 
ISM dust are clearly seen in both MiUcy Way star for- 
mation regions (top) and external galaxies (bottom). 
The top image is a color image of the Aquila star 
formation complex (Andre et al. 2010). This image 
was taken as part of the Gould Belt Survey Herschel 
Key Project, covers ~11 deg^, and was created from 
PACS 70 fim (blue), PACS 160 //m (green), and SPIRE 
250 fim (red) observations. The bottom image shows 
the complex structure of the ISM in M31 from HI ob- 
servations (green, Braun et al. 2009), embedded star 
formation using Spitzer 24 yum data (red, Gordon et al. 
2006), and unobscured young stars from the GALEX 
far-UV images (blue, Thilker et al. 2005). 



Therefore, most of the current 3D dust applications 
are intrinsically multi-wavelength in nature. 

Other difficulties are related to the complexity that 
is encountered when accessing 3D structures. The un- 
derlying grids to resolve the sink and source contri- 
butions to the radiation fields are generally discretized 
and require substantial storage and the RT calculation 
effort rises with the cell size for a given solution accu- 
racy criterion. Moreover, when modeling the structure, 
the spatial distribution of the sources and sinks has to 
be parametrized. Modeling complex structures with 
a simple spatial distribution model can lead to mis- 
leading results. Witt & Gordon (1996) showed that 
RT through a 3D fractal dust distribution is signifi- 
cantly different than through similar, but smooth dis- 
tributions. Combined with the runtimes expected for 
the 3D dust RT code, an exploration or optimization 
of the parameter space challenges the capabilities of 
current computers. Finally, more than for simple ge- 
ometries, applying 3D dust RT also struggles with the 
loss of information due to projection effects. 

As a result of non-local and non-linear effects, the 
RT equation is an integro-differential equation includ- 
ing a scattering integral; the thermal source term is 
non-linearly coupled to a double-integral equation, 
making it difficult to apply common solvers. More- 
over, the varying extinction causes changes in the 
numerical nature of the RT equation. Its character 
changes from parabolic for the diffusive transport to 
hyperbolic for freely-streaming photons; to a combi- 
nation of the two, in the numerically difficult transition 
region. Solving such a high-dimensional non-local, 
non-linear problem requires substantial computational 
resources (both computing power and memory), af- 
fecting the solution algorithms and potential limiting 
the model complexity. 

With this review we enter this evolving and dy- 
namic field to report on the significant progress that 
has been made to tackle this grand challenge prob- 
lem. Within the last 10 years, the availability of the 
high resolution images and increase computer speed 
and storage have triggered an expansion of the dust 
RT community and the development of new codes ca- 
pable of dealing with the complete 3D dust RT prob- 
lem. Many of the techniques used to solve the 3D dust 
RT problem were developed originally for ID or 2D 
geometries. The added computational complexity of 
solving the 3D problem has emphasized the need for 
highly efficient techniques, leading to refinements in 
and use of all possible 1D/2D methods in most 3D 
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codes. Applications explicitly using 3D RT include 
models of young stellar objects (Wolf, Fischer & Pfau 
1998), proto-stellar to proto-planetary disks (Indebe- 
touw et al. 2006; Niccolini & Alcolea 2006), reflec- 
tion nebulae (Witt & Gordon 1996), molecular clouds 
(Steinacker et al. 2005; Pelkonen, Juvela & Padoan 
2009), spiral galaxies (Bianchi 2008; Schechtman- 
Rook, Bershady & Wood 2012), interacting and star- 
burst galaxies (Chakrabarti et al. 2007; Hayward et al. 
2011), and active galactic nuclei (Schartmann et al. 
2008; Stalevski et al. 2012). This expansion motivated 
the "Cosmic Dust and Radiative Transfer" workshop 
held in Heidelberg in 2008 \ During the workshop, we 
realized that a review of the various techniques and ap- 
plications addressing 3D dust RT was need to commu- 
nicate common strategies between related fields. We 
felt that a review of 3D dust RT would be useful for 
coders and users of dust RT codes and people wish- 
ing to enter this field (including writers of line RT and 
(M)HD codes). 

While various solver techniques are in use for RT 
problems, 3D dust RT is commonly solved using the 
Monte Carlo (MC) technique, with some applications 
being solved using the Ray-tracing (RayT) technique. 
Since modern MC solvers make use of some RayT 
methods, this review concentrates on describing the 
spectrum of techniques based on these two approaches. 

Other RT solution methods exist but have either not 
been used to solve the complete 3D dust RT prob- 
lem in all aspects, or have shown clear disadvan- 
tages. One potential solution method is discretizing 
the RT equation, e.g. with a finite difference approach 
in spatial Cartesian coordinates and in direction, to 
create a system of linear equations (Stenholm, Sto- 
erzer & Wehrse 1991; Steinacker, Bacmann & Hen- 
ning 2002). The corresponding matrix is extremely 
difficult to solve even with powerful matrix solvers 
(van der Vorst 1992). Two more methods are used for 
RT problems, but have yet to be applied to 3D dust 
RT. The discretization can be performed on unstruc- 
tured grids (e.g., a Delaunay grid) and this can provide 
a very fast solver of RT problems. Currently, such al- 
gorithms have been updated to handle freely streaming 
photon packages and changes in the optical depth, but 
treatment of scattering has yet to be explored. Finally, 
the moment method expands the intensity as a function 
of angle using spherical harmonics as basis functions. 
It has several numerical advantages both in terms of 
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solution accuracy and storage requirements, but can 
exhibit non-physical oscillations. A common variant 
of the moment method is the Variable Eddington Ten- 
sor method (used, e.g., for 2D in the code RADICAL, 
see description in Pascucci et al. 2004). 

This review starts with the mathematical definition 
of the full 3D dust RT problem. Next, the discretiza- 
tion of the problem in spatial, direction, and wave- 
length dimensions is presented. The RayT and MC 
methods of solving the 3D RT problem are described 
in detail. Challenges in comparing RT models with ob- 
servations are discussed. A listing of existing 3D dust 
RT codes is given along with current benchmarking ef- 
forts. Finally, the review is concluded with a summary 
and discussion of the future of 3D RT. 

2. THE 3D DUST RADIATIVE TRANSFER 
PROBLEM 

2.1. The radiative transfer equation 

The stationary radiation field is described by the 
specific intensity I(x, n, A), where x gives the location 
in space, « is a unit vector indicating the direction of 
the radiation, and A its wavelength. The specific in- 
tensity represents the amount of energy carried by ra- 
diation in a unit wavelength interval, which is trans- 
ported per unit solid angle and per unit time across an 
element of unit area perpendicular to n.* The contin- 
uum radiative transfer equation (RTF) describes how 
the specific intensity varies as a result of interactions 
with a medium filled with sources and sinks. In its gen- 
eral form, it can be written as (see e.g Chandrasekhar 
1960; Rybicki & Lightman 1979) 

n VIix, n. A) = -k(x. A) p(x) I(x, n, A)+j{x, n. A). (1) 

The left-hand side of this equation represents the 
change of the intensity over an infinitesimal distance 
along the path determined by the position x and the 
propagation direction n. The first term on the right- 
hand side represents the extinction, i.e. the loss of 
radiant energy, when radiation passes through matter 
Here k{x. A) is the mass extinction coefficient, and p(x) 
is the mass density. The second term on the right-hand 
side represents the source term, i.e. the new luminosity 



The specific intensity can be defined as tlie intensity per unit of 
wavelengtli or per unit of frequency. The convention chosen is typ- 
ically indicated by a subscript, i.e. as /,i or ly. We adopt the per 
unit of wavelength convention, and drop the subscript in order not to 
overload the notations. 
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released into the medium at x in the direction n. The 
complexity of the RTE depends on the nature of the 
source and sink terms, i.e. the different physical pro- 
cesses that are responsible for extinction and emission. 

An alternative form of the RTE (1) explicitly uses 
the distance s along the path defined by a position x 
and propagation direction « as a variable. We then ob- 
tain 



dl 



(s. A) = -k(s, A)p{s) I(s, A) + j(s, A). (2) 



If we assume for now that the source term j does not 
depend on the intensity /, we can readily solve this 
differential equation, 



I(s,A)^ f j(s',A)e-'''-'''''-^'>ds' , 

\J —OQ 



(3) 



with the optical depth between two positions defined 
as ^ 

t{si,S2,A)^ I K{s,A)p{s)ds. (4) 

J Si 

Expression (3) has a simple physical interpretation: it 
expresses that the intensity at any position s along a 
path results from the emission at all anterior points s' 
along the path, reduced by a factor e"'"*'"' ''*' to account 
for the extinction by the intervening matter. It is impor- 
tant to stress that the expression (3) is only a. formal so- 
lution of the RTE, and not a very useful solution of the 
RTE. Indeed, the emissivity generally does not only 
depend on position, dkection and wavelength, but also 
on the specific intensity itself. The formal solution (3) 
is then no more than an integral equation version of the 
radiative transfer equation itself; this is particularly the 
case for dust radiative transfer In the following sub- 
sections, we gradually build the dust RTE by including 
the various relevant physical processes. 

2.2. Primary emission and absorption 

Two important and obvious processes to take into 
account in a dusty medium are primary emission and 
absorption. Primary emission accounts for the radia- 
tive energy added to the radiation field - this is often 
stellar emission, but it can include, for example, radia- 
tion from an active galactic nucleus, emission line ra- 
diation from ionized gas or Bremsstrahlung. In general 
form, it can be characterized by a function n, A). 
Absorption is the process in which electromagnetic ra- 
diation is taken up by dust grains and transformed into 
the internal energy. It is characterized by the absorp- 
tion coefficient k^\,^\ for a given chemical composition. 
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size and shape of a dust grain, the absorption coeffi- 
cient can in principle be determined at any wavelength 
(Mie 1908; Purcell & Pennypacker 1973; Draine 1988; 
Mishchenko, Travis & Lacis 2002; Min, Hovenier & 
de Koter 2005). 

When we take only primary emission and absorp- 
tion by dust into account, the RTE (1) becomes 

d/ 

-{x,n,A) - -K2b^{x,A)p{x)I{x,n,A) + jJx,n,A). 

(5) 

This equation is a simple first-order differential equa- 
tion and the solution can be found by just integrating 
along the line of sight, as for the formal solution (3). 
For general 3D geometries, this integration is done nu- 
merically. 

2.3. Including scattering 

The complexity of the RTE increases substantially 
when we take scattering into account. Scattering, as 
absorption, removes radiation from a beam and hence 
accounts for a second sink term in the RTE, the ef- 
ficiency of which is quantified by the scattering co- 
efficient /Csca- Rather than converting the radiation to 
internal energy, it emits the same radiation in a dif- 
ferent direction. Scattering hence does not only im- 
ply a second sink term, but also a second source term. 
The scattering phase function <i>{n, n' ,x. A) describes 
the probability that a photon originally propagating in 
the direction n' and scattered at the position x, will 
have n as its new propagation direction after the scat- 
tering event. Given this definition, the phase function 
satisfies the normalization 

I <i{n,n',x,A)dQ.' ^ I <i)(n,n',x,A)dn ^ I. 

JAn JAti 

(6) 

Adding the sink and source terms due to scattering to 
the RTE, we obtain 



d/ 

As 



{x, n. A) - -Kf.^t(x, A) p(x) lix, n. A) + jt{x, n. A) 



+ Kscii(x,A)p(x) I <i>(n,n' ,x,A) I(x,n' ,A)dD.' . 



(7) 



where the extinction coefficient /Cext — ''^abs + ""sea- 

Un- 
like the simple differential equation (5), this equation 
is an integrodifferential equation in which the radia- 
tion fields at all different positions and in all different 
directions are coupled. 
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In many RT calculations, scattering is important 
and adds significant complexity given that dust scat- 
tering is anisotropic. This is especially true for UV 
to NIR wavelengths: observations indicate that the 
dust scattering albedo at these wavelengths is at least 
50%, and that scattering off dust grains is strongly 
anisotropic (Witt et al. 1992; Calzetti et al. 1995; 
Burgh, McCandliss & Feldman 2002; Gordon 2004)^ 
Using an isotropic phase function or other approxima- 
tions such as an isotropic two-stream approximation 
or an (effective) forward scattering (e.g. Code 1973; 
Natta & Panagia 1984; Calzetti, Kinney & Storchi- 
Bergmann 1994) might not always be physically jus- 
tifiable. Several authors demonstrated that an im- 
proper treatment of anisotropic scattering leads to sig- 
nificant errors (e.g. Bruzual, Magris & Calvet 1988; 
Witt, Thronson & Capuano 1992; Baes & Dejonghe 
2001). Even for radiation at MIR wavelengths or 
longer, where the scattering off common ISM grains is 
low, an application requiring the heating of the grains 
to be determined will need to properly handle scatter- 
ing (Nielbock et al. 2012). 

The Henyey-Greenstein (HG) phase function (Henyey 
& Greenstein 1941) is the most widely used parametriza- 
tion of the dust phase function and provides a good 
single parameter approximation. Dust grain models 
predict small deviations from a HG phase function and 
other parametrizations or numerical phase functions 
can be used for increased accuracy (Kattawar 1975; 
Hong 1985; Draine 2003b). 

2.4. Radiative transfer in dust mixtures 

In any realistic dust medium, there is a range of dif- 
ferent types of dust grains present, with various chem- 
ical compositions, sizes, shapes and number densi- 
ties. Each of these different grain types / is charac- 
terized by its own absorption coefficient K-^bsA^), scat- 
tering coefficient /fsca,i(^) and scattering phase function 
<!),■(«,«', /i). If we denote the relative contribution of 
each grain type / at the position x to the total dust num- 



^ Up-to-date versions of the Gordon (2004) plots of albedo 
and scattering phase function asymmetry versus wavelength can 
be found at http://www.stsci.edu/-kgordon/Dust/Scat_ 
Parain/scat_data . html. 



ber density as Wi{x), the transfer equation becomes 

^{x, n,A) = - V Wi(x) K^^tjiA) P(x) I(x, n, A)+jAx, n. A) 
as ^ 
I 

+^ Wi{x) K,^^,i{A) Pix) I 1>,(«, A) I(x, n'. A) AO! . 

: Jin 



(8) 



It is straightforward to see that (8) is formally identical 
to equation (7) if we define the absorption coefficient, 
scattering coefficient, extinction coefficient and scat- 
tering phase function of a dust mixture as 

K^uix, /I) = ^ Wi(x) iW, (9a) 

/ 

'''sea, iiX), (9b) 

/ 

K^^liX, /I) = ^ Wi(x) i(X), (9c) 
i 

^, , Y.i^iix)K,^^,i{A)^iin,n' ,A) 

0(/i, n , jc, /I) = — — . (9d) 

Li Wi(x) K^ci^M) 

As long as we only consider primary emission, ab- 
sorption and scattering, radiative transfer in dust mix- 
tures is hence completely identical with radiative trans- 
fer in a dust medium with a single average dust grain 
and no approximations are needed (Martin 1978; Wolf 
2003a). 

2.5. Including dust emission 

Apart from primary emission, absorption and scat- 
tering, a fourth physical process to be taken into ac- 
count in dust RT is the thermal emission by the dust 
itself. Dust grains that absorb radiation re-emit the ac- 
quired radiative energy at wavelengths longwards of 
about 1 yum. To account for this astrophysical process, 
we need to incorporate the third source term jdix. A) in 
our RTE, 



d/. 
As 

+ Kse^{x,A)p{x) 



(x, n, A) - -Kf-^tix, A) p(x) I(x, n, A)+ jt(x, n, A)+ ji(x, A) 

I ^(n, n',x. A) I(x,n', A) dQ.' . 

(10) 



It might seem as if the dust emissivity term is a simple 
extra source term similar to the primary stellar emis- 
sivity term. Its exact form is strongly dependent on 
which physical emission processes are important, and 
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often the dust emissivity term depends in a compli- 
cated and non-linear way on the intensity of the ra- 
diation field itself 

A common assumption is that the dust grains are in 
thermal equilibrium with the local interstellar radiation 
field. In this case, the emissivity of the population of 
grains of type / can be written as a modified blackbody 
emission characterized by an equilibrium temperature 
Ti{x). Summing over all populations, we obtain 

Mx, ^) = 2 Wi(x) K,ksAA)p{x) B[Ti(x), A), (11) 

with B{T, A) being the Planck function. The equilib- 
rium temperature of each type of grain is determined 
by the balance equation, i.e. the condition that the total 
amount of energy absorbed equals the total amount of 
emitted energy, 

XOO y-»00 
K.h.AA) J(x, A)dA^ J ^,bs,,(^) B^{Ti(x), a) dA, 

(12) 

where J{x, A) represents the mean intensity of the ra- 
diation field. 



J(x, A) 



= — r lix, n. A) dn. 

4^ J4n 



(13) 



It is important to note that the equilibrium tempera- 
ture of the dust grains depends explicitly on their size 
and chemical composition. At the very same loca- 
tion, dust grains with different sizes and/or chemical 
compositions will obtain different equilibrium temper- 
atures. So far, we have easily combined the absorp- 
tion and scattering due to different kinds/sizes of dust 
grains in the RTE without any approximations. The 
same cannot be done for the thermal re-emission term. 
One might be inspired by equations (9) and average the 
temperatures of the different grains to a "mean" tem- 
perature. This results in reducing the complexity of 
the dust mixture at a given position to a single mean 
grain that will reach a single equilibrium temperature. 
While this could be useful, sufficient or even necessary 
for certain applications, it is important to be aware that 
this is a simplification of the RT problem that is not 
physically correct (e.g. Wolf 2003a). 

While the assumption of thermal equilibrium is use- 
ful in some applications, it breaks down in others. Par- 
ticularly important is the case where the dust medium 
contains very small dust grains (including polycyclic 
aromatic hydrocarbons [PAHs]). Large dust grains 
reach thermal equilibrium and emit as modified black- 
bodies with an equilibrium temperature. However, 



small dust grains have small heat capacities, and the 
absorption of even a single UV/optical photon can sub- 
stantially heat the grain. These small grains will not 
reach an equilibrium temperature but will instead un- 
dergo temperature fluctuations that lead to grain emis- 
sion at temperatures well in excess of the equilibrium 
temperature. The emission from small grains is nec- 
essary to explain the observed mid-infrared emission 
of many objects (e.g. Sellgren 1984; Boulanger & Per- 
ault 1988; Helou et al. 2000; Smith et al. 2007; Draine 
et al. 2007). When including such transiently heated 
dust grains, the dust emissivity changes from expres- 
sion (11) to 

r y-»tX) 

Mx, A) = y wi(x) K,^,,/A)p(x) P;{T, X) B(T, A) dT 

(14) 

Here Pi{T,x) is the temperature distribution for dust 
grains of type / at the location x. This temperature 
distribution depends on the chemical composition and 
size of the dust grains, as well as on the intensity and 
hardness of the radiation field in which it is embedded. 
Several methods have been developed to calculate the 
temperature distribution of small dust grains, using ei- 
ther matrix operations or time averages (Dwek 1986; 
Desert, Boulanger & Shore 1986; Guhathakurta & 
Draine 1989; Siebenmorgen, Kruegel & Mathis 1992; 
Draine & Li 2001; Compiegne et al. 201 1). The result 
is that the dust source term is an intricate, non-linear 
function of the specific intensity, which adds a signifi- 
cant complexity to the RT problem. 

Finally, thermal emission is not the only emis- 
sion process of dust grains: additional non-thermal 
processes are extended red emission (Witt & Boro- 
son 1990; Smith & Witt 2002) and blue luminescence 
(Vijh, Witt & Gordon 2004). These non-thermal pro- 
cesses can account for a substantial fraction of the sur- 
face brightness of interstellar clouds at optical wave- 
lengths (e.g. Gordon, Witt & Friedmann 1998; Witt 
et al. 2008). Both processes can be included as an 
additional term in the dust source term j^ix. A) in the 
RTE (10). 

2.6. Radiative transfer of polarized radiation 

The specific intensity I(x, n. A) is not a complete de- 
scription of the radiation field, as it only describes un- 
polarized light. Scattering of photons from dust grains 
naturally produces polarized radiation (e.g.. Fig. 2 of 
Gordon et al. 2001). In addition, aligned dust grains 
also produce polarized radiation (Whitney & Wolff 
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2002). While alignment of grains has been demon- 
strated observationally for many decades, the physical 
mechanism for grain alignment is still matter of debate 
(Lazarian 2007). 

The most common way to describe polarized radi- 
ation makes use of the Stokes vector S = (/, Q, U, V), 
where / is the total specific intensity, Q and U the lin- 
early polarized intensity in two axes rotated by 45 de- 
grees one from the other, and V the circularly polarized 
intensity. The four components do not form a preferred 
basis of this space, but were chosen because they can 
easily be measured or calculated. 

The RT formalism we have described in the previ- 
ous subsections can be extended to include polarized 
radiation. Instead of a single RTE, we then obtain 
a vector RTE, or equivalently, a set of four coupled 
scalar equations. 



d5 



(x, n. A) = -K^^iix, A)p{x) S{x, n, A)+j^(x, n, A)+j^(x, A) 



+ K^c^{x,A)p(x) I M(n,n',x,A)S(x,n',A)dD.' . 

(15) 

A first complication is the scattering source term, 
where the phase function <!>(«, n' ,x. A) is now replaced 
by a 4 X 4 scattering (or Mueller) matrix M{n, n',x. A) 
which describes the changes in the Stokes vector when 
radiation is scattered from propagation direction n' to 
a new propagation direction n. For a full description 
see e.g. Bohren & Huffman (1983), Fischer, Henning 
& Yorke (1994) or Code & Whitney (1995). When 
we consider RT of polarized light through a dust mix- 
ture, each type / of grains is characterized by its own 
Mueller matrix Aii{n, n' , A). The radiative transfer 
equation can then still be written as (15), with 

, I,iWiix)K,c.^jiA)Miin,n',A) 
M(n,n,x,A)= 7- ttt ■ (16) 

A second complication is the thermal emission term: 
for thermal emission from aligned grains the full 
Stokes vector is used in the emission. 

3. THE DISCRETE 3D DUST RADIATIVE TRANS- 
FER PROBLEM 

A general analytical solution of the stationary 3D 
dust RT equation is not possible for any of the non- 
symmetric applications mentioned in the introduction. 
To apply numerical solution techniques, solvers gener- 



ally discretize the solution vector or the physical prop- 
erties in the RTE. The quantities requiring discretiza- 
tion are the three spatial coordinates, the two direc- 
tional coordinates, the wavelengths, and/or the dust 
properties. 

A major concern when solving a 6D integro- 
differential equation is the size of the solution vector 
With a resolution of 100 points in each variable, the 
intensity vector has 10'^ entries. Handling this inten- 
sity vector requires an enormous amount of computer 
memory and speed. Currently, many solution algo- 
rithms avoid this requirement by not storing the full 
solution vector. There are applications where the full 
direction-dependent intensity is needed (e.g., when the 
radiation pressure impacts the gas kinematics). Due to 
this high dimensionality of 3D RT, the choice of appro- 
priate grids is mandatory to effectively apply existing 
solution techniques and minimize memory usage and 
runtime. The solution techniques used does influence 
the choice of grids (e.g., RayT versus MC). 

Another concern is that for a discrete solution vec- 
tor, the physics is only solved at the grid resolution 
even if the solver used has no intrinsic error. Thus, 
if the grid is too coarse, a strong change in the inten- 
sity due a change from optical thin to thick may not 
be resolved and the derived intensity would have large 
systematic errors. This makes the choice of the grid 
crucial. 

3.1. Spatial grids 

For many astrophysical applications, the density 
values cover orders of magnitudes and are highly 
structured (e.g., a turbulent gas cloud with filaments 
and shocks). The same is true for the sources which 
can be dust grain emission, the distribution of stars, 
or a layer within a photon dissociation region (PDR) 
emitting in the MIR. Consequently, the complexity 
of spatial grids in 3D continuum RT ranges from 
simple linear Cartesian grids (Stenholm, Stoerzer & 
Wehrse 1991) to adaptively refined Cartesian grids 
(Kurosawa & HiUier 2001; Niccohni & Alcolea 2006; 
Lunttila & Juvela 2012) to multi-wavelength AMR 
grids (Steinacker, Bacmann & Henning 2002). Com- 
plex dust distributions are illustrated for three different 
cases in Fig. 2 using refined Cartesian grids. In princi- 
ple, the RTE could be solved on unstructured, dynamic 
grids like those used in line RT (Petkova & Springel 
2011; Paardekooper, Kruip & Icke 2010). Finally, the 
density or source distribution could be given using an 
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Fig. 2. — Examples of advanced spatial grids that are 
currently being used for 3D dust RT calculations, from 
Saftly et al. (2012), are shown. The examples are 2D 
plane cuts through octtree-based grids used in Monte 
Carlo simulations; similar grids are applied in ray- 
tracing techniques. On each row, the left panel shows 
the dust density, the central panel the grid distribu- 
tion in a regular octtree structure, and the right panel 
the grid distribution in a barycentric octtree structure. 
The top row represents a regular, analytical model of 
a double-exponential disc with a three-armed logarith- 
mic spiral density perturbation. The middle row corre- 
sponds to a clumpy AGN model, consisting of a large 
number of optically thick clumps in a smooth density 
distribution (Stalevski et al. 2012). The bottom row 
corresponds to a late-type disc galaxy model resulting 
from N-body/SPH simulations with strong supernova 
feedback (Rahimi & Kawata 2012). In all cases, the 
grids contain between 3 and 4 million cells. 



analytical formula and the 3D dust RT would be a 
complex function with many parameters. 

The optimal grid would be based on changes in the 
radiation field intensity. As we don't know the radi- 
ation field a priori (this is the goal of the RT calcula- 
tions), it is extremely challenging to generate such an 
optimal grid. First, the spatial grid is only 3D while 
the intensity is defined in 6D (spatial, direction, and 
wavelength). For example, the intensity in a certain 
direction can remain constant between two positions, 
while the intensity in another direction can vary dras- 
tically between the same positions. The radiation field 
also depends on the wavelength, so the optimal spatial 
grid is different for each wavelength. The dust mass 
density, or the optical depth, could serve as the starting 
point on which the grid could be based (see e.g. Fig- 
ure 3, which shows 3D RT octree grids with a refine- 
ment criterion based on the total dust mass in a cell). 
However, the cell optical depth alone is not suitable 
to define a grid refinement criterion, as the strength of 
the radiation field can show strong gradients even in 
regions where the optical depth is small. In summary, 
the optimal grid should combine the details of the dis- 
tribution of both the dust and source terms such that it 
captures the variation of the radiation field intensity. 

As a result of the difficulty in creating an optimal 
spatial grid, various approximations are used. To re- 
view the variety of grids and their purpose, it is con- 
structive to distinguish three spatial grid classes. They 
are used in various combinations by different solvers 
and in the different astrophysical communities. 

Local mean intensity storage grids 

This most common class of grids stores the mean 
intensity with the goal to resolving its variation. In the 
end, the grid's resolution determines the spatial reso- 
lution of the obtained solution. 

To achieve a reasonable sampling of the gradients 
in the mean intensity, it is possible to calculate a series 
of spatial grids that are refined using the local opti- 
cal depth averaged over all directions for each wave- 
length. The computational effort to calculate such 
grids is negligible compared to the effort of solving the 
RT equation (Steinacker, Bacmann & Henning 2002). 
The drawback of multi-wavelength grids is that a large 
number of interpolations have to be performed in or- 
der to assemble the wavelength dependent radiation 
field. Such interpolations are time consuming and in- 
troduce interpolation errors. As a compromise, grids 
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optimized for tracking the variation in the local optical 
depth for the wavelengths that dominate the radiative 
energy locations have the advantage of resolving the 
important regions at the expense of excess grid points. 
An alternative way of building the grid was proposed 
by Niccolini & Alcolea (2006) where an initial calcu- 
lation is used to define the explicit location of a num- 
ber of photons and these locations then used to refine 
the grid to have a constant absorption rate in each grid 
ceU. 

Density and source grids 

The grids storing density and source information 
can originate from a discretization of simple physical 
models to deep structure grids designed to resolve the 
kinematic processes. The latter one is often realized 
as adaptive mesh refinement grid or smoothed particle 
(SP) distributions used in hydrodynamical simulations 
(e.g. Steinacker et al. 2004), typically incorporating 
several tens of millions of cells and many levels of re- 
finement or SPH particles. While RT calculations can 
be performed directly using such density grid (Abel & 
Wandelt 2002), in most cases the density information 
is stored on coarser grids to meet the storage and speed 
limitations of the 6D RT. 

Solution grids 

The third class of grids is designed to calculate 
the solution directly at the grid points by advancing 
from one cell to the next in one step. The refine- 
ment criterion is defined to minimize solution errors 
of a certain order (Steinacker et al. 2002) or to provide 
a flexible grid that minimizes solution errors due to 
the coarse spatial coverage of the physical quantities 
(Paardekooper, Kruip & Icke 2010). These grids are 
well suited for (M)HD codes. Solving the RT equa- 
tion directly on a grid accumulates discretization er- 
rors causing, e.g., a smearing of beams in the case of 
finite differencing solvers. Additional numerical meth- 
ods or higher-order corrections are used to compensate 
for these numerical diffusion errors. It should be noted 
that in some applications the source function may also 
contain a smaller number of discrete sources such as 
stars. These can be considered outside the main spatial 
grid or using smaller subgrids (Stamatellos & Whit- 
worth 2005). 



3.2. Direction grid 

There are two major challenges in defining a fixed 
direction grid for 3D RT. First, the radiation field can 
be strongly peaked due to nearby sources, and a smear- 
ing of this beam due to insufficient direction resolu- 
tion will result in remote regions not being illuminated 
accurately. Since the optimal intensity discretization 
can vary greatly across the domain, a globally opti- 
mal refinement of direction space is not usually pos- 
sible. Second, even direction grids optimized to be 
equally-spaced on the unit sphere (Steinacker, Thamm 
& Maier 1996) require a large number of points pro- 
vide good angular resolution. Another possible reg- 
ular direction grid is given by the HEALPix method 
which subdivides the unit sphere in pixels of equal 
size (Gorski et al. 2005). For example, a grid with 
600 direction points equally spaced on the unit sphere 
provides a mean resolution of ~10° only, and it takes 
10000 grid points to achieve a mean resolution of 
~2.5°. For a protoplanetary disk, this resolution cor- 
responds to a hot dust region with the size of 4!AU 
placed at a distance of 100 AU from the central star. 

It is important to note that the two main RT solution 
techniques treat the discretization of directions quite 
differently. In MC, the direction of the photons is not 
discretized, but sampled from a probability function. 
For example, in MC the calculation of the scattered in- 
tensity is based directly on the scattering phase func- 
tion, allowing an arbitrarily precise solution. In RayT, 
a discrete direction grid is used for all calculations. 
The scattered intensity calculation is done on the di- 
rection grid, potentially undersampling the scattering 
process in the direction space of the solution. 

Another issue in direction space is the divergence of 
rays or photon tracks. The chance to miss physically 
important parts of the domain increases with the dis- 
tance from the current point or cell. When tracing the 
radiation from a single source, this can lead to large er- 
rors in the computed radiation field for distant cells, or 
long runtimes to increase resolution using more pho- 
tons or rays. Solutions to this issue exist and are dis- 
cussed in the MC and RayT sections. 

3.3. Wavelength and dust grain grids 

It is important to consider variations in source spec- 
tra and dust opacity when choosing a wavelength grid. 
The spectrum of the radiation sources should covered 
well enough to resolve the overall shape and any im- 
portant spectral features (e.g., emission lines). The 
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wavelength grid should resolve variations in the grain 
properties (e.g., opacities and scattering properties). 
Where there are features in the dust properties (e.g., 
2175 A bump, MIR aromatic/PAH features), the wave- 
length grid should resolve the feature, ideally includ- 
ing a point at the maximum of each feature as well as 
enough points to Nyquist sample the feature. 

Beside the discretization of the variables entering 
the intensity directly, the size distribution of the grains 
needs to be discretized if not given analytically. The 
grain size discretization can have a strong impact on 
the radiation field. The extinction of the radiation is 
the sum of the extinction of the different species, but 
the emissivity of individual grains is coupled to the in- 
tensity of the incoming radiation field (see Sec. 2.3). 

4. THE RAY-TRACING SOLUTION METHOD 

Ray-tracing (RayT) is a method widely applied in 
physics and computer graphics to describe the propa- 
gation of electromagnetic waves or particles through 
a medium with varying properties. Important appli- 
cations outside astrophysics cover the propagation of 
radio signals in the ionosphere, the investigation of 
heating by plasma waves, sound waves in the ocean, 
optical design of lenses, tomographic reconstruction 
of the Earth's interior, and image generation in com- 
puter graphics. In radiative transfer, RayT is used to 
follow the change of intensity in a particular direc- 
tion which is the basic transport problem described by 
Eq. (1). The MC solution technique described in the 
next section is a sophisticated variant where a prob- 
abilistic approach is taken to choose the direction of 
the photon package representing the ray. RayT solvers 
have been used for a number of 2D continuum RT ap- 
plications (see e.g. the benchmark comparison in Pinte 
et al. 2009) and 3D RayT is based on many techniques 
first developed in 2D. 

In the following subsections, we describe the ba- 
sic ingredients and challenges when using RayT as a 
solver of the general 3D continuum RTE. The solution 
on a single ray under various numerical and physical 
conditions is given in §4.1. The global solution and 
the treatment of source terms coupling directions are 
discussed in §4.2. 

4.1. RayT Solution for Single Ray 

The elementary RayT operation is to solve the first- 
order differential equation (2) within a spatial density 
grid cell along a given direction. The mass density 



po, the mass extinction coeflicient kq, and the source 
function jo are assumed to be constant within the cell 
for a given wavelength, so that we can determine the 
intensity I{s + As, A) from I(s, A) using Eq. (3) 

I(s + As, A) = lis, A)e-^<'^'^ + ^^^^ (l - e-'«^'>) 

To(A) ^ ' 

(17) 

with To(/i) = Kf){X)p{^As. For a ray crossing several 
cells, the numerical task is to determine the cells that 
are hit by the ray, calculate the intersection points 
with the cell borders, and use Eq. (17) to calculate the 
change in intensity along the ray in each cell. The first 
step can be time-consuming, since an adaptively re- 
fined grid is often stored as a tree, and neighbor-search 
calculations are required. The error in the intensity is 
solely defined by the the finite resolution of the under- 
lying spatial grid. 

For clarity in the notations used, we note that the 
solution of the RTE along rays with constant direction 
is also called the method of short or long character- 
istics. Long characteristics refers to calculating the 
radiation field along a ray through the entire compu- 
tational domain. Since several rays can cross a certain 
cell causing redundant calculations, pre-calculated lo- 
cal column density or optical depth values can be used 
(short characteristics) to interpolate the contributions 
along a ray. Usually, a sweep of the ordering of the 
grid is required to ensure that, for a given direction, all 
information about positions before the currently con- 
sidered point are given before performing the step de- 
scribed by Eq. (17). The method is less accurate than 
long characteristics due to the accumulation of inter- 
polation errors. For combined applications of RT and 
(M)HD, hybrid methods have been proposed combin- 
ing the advantages of short and long characteristics 
(Rijkhorst et al. 2006). 

Beyond the Spatial Grid Resolution 

There are cases where the best RayT solution is not 
performed at the density grid resolution. First, the spa- 
tial resolution of the density from an Adaptive Mesh 
Refinement (AMR) MHD calculation with many re- 
finement levels can be too fine for a RT calculation 
to be done in a reasonable time. One way to deal 
with this case is to interpolate the density to a coarser 
grid, but this loses some of the information on the 
physical structure obtained in the AMR run. Second, 
the density grid shows strong gradients. One way to 
soften the gradients is to interpolate the density to a 
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finer grid. But this can still leave the density changes 
abrupt. Instead an interpolation scheme can be used 
to find the density in between grid cells. Finally, the 
density is given analytically and no step size infor- 
mation is a available. Note that in practice analytical 
density structures are not automatically simpler than 
AMR density grids. The density used for the RT mod- 
eling of Spitzer images of the molecular cloud LI 83 by 
Steinacker et al. (2010), e.g., involved 100 3D clumps 
with Gaussian profiles and 700 free parameters. 

The simplest approach to solve the RTE in these 
cases is to apply an upwinding first-order finite dif- 
ference scheme. The relation between two intensities 
located at s and 5 -H Ai along the ray is then 

I(s + As, A) = I{s, A)[l- t{s, As, A)] + j(s, A)As. (18) 

using the local optical depth t{s, As, A) = k(s, A)p(s)As. 
The step size As is chosen to be small enough to stay in 
the optically thin limit allowing the Taylor expansion 
of the exponential to be truncated after the first term 
(exp[-T] ~ 1 -t). The advantage of this scheme is that 
it is fast, the disadvantage is that first-order errors can 
accumulate along the ray. Moreover, the steps become 
very small in regions of high optical depth, although 
the radiation field takes a simple form in this limit. 

The concept of "trial" steps can improve the accu- 
racy. Using a Runge-Kutta scheme, a step with the size 
As/2 can be made to calculate a second estimate of the 
derivative. The scheme then advances to the next point 
by using a linear combination of both derivatives. The 
linear factors are chosen by comparison with the Tay- 
lor series of I{s) and letting the first-order term vanish. 

This improvement can be repeated to achieve a bet- 
ter solution, at the cost of repeated calculations or 
look-up of the density and source terms. An advanced 
ray-tracer based on a 5th-order Runge-Kutta accuracy 
has been proposed by Press et al. (2002). It is coupled 
with an adaptive step size control using an embedded 
form of the 4th-order Runge-Kutta formula. As pa- 
rameters for the error truncation, values determined by 
Cash & Karp (1990) are used. The tracer steps are cho- 
sen adaptively and with full error control. Such a tracer 
is therefore the first choice for astrophysical RT prob- 
lems with moderate optical depth variations which are 
not solved at the resolution of the density grid. This 
explicit error control is an important characteristic of 
the RayT solution technique. 



High optical depths 

There are 3D RT applications with strong gradients 
in the source function or in the local optical depth. In 
particular, the case of high optical depth r » 1 occurs 
in all star formation regions as well as in AGN tori. All 
algorithms used to solve the RT encounter problems in 
correctly describing the intensity changes in the opti- 
cally thick case. Adaptive grid RT solution techniques 
refine the regions into too many cells. Only moment 
methods which are applied in conjunction with (M)HD 
solvers are well-posed to treat the optical thick regime, 
but in turn fail to describe a strongly peaked radiation 
field arising in low optical depth regions. 

To illustrate the difficulties of high optical depths 
for the RayT method, we assume a simple modi- 
fied black body thermal source term for a single dust 
species j(s,A) = k(s,A)p(s)B{s,A) in Eq. (18), which 
becomes 

I(s + As, A) - I{s, A) = t(s, A) [B{s, A) - I(s, A)] . (19) 

When the ray moves from one to the next point, most 
of the radiation at s is damped, as is the source con- 
tribution along the path, so that the local radiation at 
s -H As is dominated by the source contribution. Hence 
B - I is small, suggesting that the solver can perform 
large steps. But B - / is multiplied by the large opti- 
cal depth T amplifying any change in the source term 
that arises from the spatial variation in the temperature. 
Thus, to meet accuracy requirements, the tracer must 
perform small steps. 

The second reason to modify the ray-tracer is the 
two exponential functions entering the radiation equa- 
tion. First, the exponential containing the optical depth 
has to be calculated precisely along the ray. Given a 
limited computational range of the computer, an op- 
tical depth of 1000 usually exceeds these limits and 
makes it necessary to renormalize the intensity. Sec- 
ond, the Planck function rises sharply with A for tem- 
peratures in the Wien part T(x) < hc/Ak. 

To solve the two exponential problem, Steinacker, 
Bacmann & Henning (2006) proposed to use the trans- 
formation to the relative intensity D - 1/(1 + B). For a 
vanishing source function, it approaches unity, and for 
the optically thick part with I{s, A) ^ B(s, A), D(s, A) 
has a value close to 1 /2. The authors showed that the 
transformation avoids numerical problems caused by 
the exponentials and that the intensity error amplifica- 
tion by the transformation is less than a factor of 2. As 
criterion for the use of an approximate solver in the 
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optically thick region D(s, A) ^ I /2, they derived 
\B(s2,A)-B(suA)\ 1 



T(Si,A) 



< e 



(20) 



where e is a positive limit for the relative difference 
\B - I\/B. The condition is fulfilled either due to small 
relative changes in the source term or a large local op- 
tical depth. In a pre-calculation along the ray, the re- 
gions where the condition is fulfilled can be identified. 
Then, a fast solution is obtained for these region by 
performing large steps while assuming D(A, x) - I /2. 
Applying the scheme to the case of a massive molecu- 
lar cloud core illuminated by a nearby star, Steinacker, 
Bacmann & Henning (2006) verified speed-up factors 
of several hundreds compared to fourth-order Runge- 
Kutta solvers. 

4.2. Ray location and RTE global solution 

After discussing the solution methods for a single 
ray through the computational domain, the next step is 
to determine how to place the rays to ensure that the 
radiation field is correctly calculated. This is a critical 
part of the solution process. 

Thermal Emission 

We describe how to place the rays in order to cal- 
culate the intensity of the radiation field when the 
source term is dominated by thermal emission from 
dust grains. The various ray patterns used in RayT 
are illustrated in Fig. 3. For what follows, the term 
"placing a ray" will be used for the basic RayT step 
that includes defining the ray by a point in space and 
a direction, solving the intensity along the ray as pre- 
viously described, and storing the absorbed energy in 
the cells that are crossed. According to Eq. (17), the 
intensity of the absorbed radiation per cell is 

lis. As, A) = lis, A)(l- e-''>^'A+MA)As 1 - ^ (l ■ 

(21) 

Its contribution to the local radiation field then is cal- 
culated from the formula for the mean intensity (13) 
and the balance equation (12). In each RayT step, the 
local source term (11) contains the thermal energy of 
the currently crossed cell; this is updated with each ray 
crossing. 

In RayT, pre-calculation steps are done to analyze 
the specific RT problem and accelerate the computa- 
tions. The optical depth affects the transport within 



the domain, the thermalization of the radiation, and the 
appearance of the object on the t(A) ^ 1 -surface for 
the observer. Therefore a penetration depth analysis is 
performed at all wavelengths determining the optical 
depth distribution with respect to all discrete sources as 
well as to the observer (see, e.g., Steinacker, Bacmann 
& Henning 2002). For this calculation, rays are placed 
from the sources to all grid cells (Fig. 3a,b). In addi- 
tion, we calculate the optical depth from the source to 
the cell using 

Tsouisi ,S2,A)^ I k{s, A)pis)ds = ^ Ki(A)piAsi 

(22) 

with the ray crossing A^^ cells with their individual /c,-. 
Pi, and crossing length A,?,. Tsou ~ 1 marks the region 
where most of the source radiation is reprocessed. 

The solver also calculates the optical depth from 
each cell to the observer Tobs as described above 
(Fig. 3f). This information is used to understand which 
regions are shielded at which wavelengths from the ob- 
server, and to calculate the final images once the main 
RayT has been performed. 

The second pre-calculation is to propagate the ini- 
tially deposited discrete source energy and the radia- 
tion field at the domain boundary through the domain. 
The calculation is done on a regular grid of rays like 
the one shown in Fig. 3e. This calculation also deter- 
mines the optical depth between cells, t^c, on a coarse 
spatial grid at all wavelengths, in this way providing 
information on which regions in the domain exchange 
significant amounts of radiation. In star formation ap- 
plications, e.g., regions often do not connect signifi- 
cantly at UV wavelengths. 

If the source of radiation is the interstellar radiation 
field, the pre-calculation can be very time-consuming, 
as it comprises rays from all directions to all cells 
a.t all wavelengths for which the optical depth can 
-To(/i)y;ach 1. However, there are many cases where this 
pre-calculation is less time-consuming. For example, 
when the dust properties are constant in the domain as 
then the optical depth is t^ou = kNcoI and this is an 
integral only over the density. 

The main RayT calculations are performed sepa- 
rately for each wavelength. Rays are placed through 
the grid points in the Tsou ~ 1 -layers covering all 
directions. In this way additional resolution is pro- 
vided where the largest changes in the radiation field 
are expected. Moreover, additional rays from discrete 
sources are placed to help resolving the illumination 
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of the Tsou ~ 1 layers. This can be important, e.g., 
in the case of an accretion disk atmosphere that is il- 
luminated by a protostar through a narrow solid an- 
gle. It is important to note that the order in which 
the computation is done for each wavelength influ- 
ences the convergence. In RayT, information has al- 
ready been transported from the sources to the cells 
during the pre-calculation step. Therefore, starting 
with wavelengths covering the peak of the re-emitted 
photon energy (e.g., in the MIR) can speed up the in- 
formation transport in thermally dominated problems. 
For problems dominated by low optical depths in the 
UV/optical, the order is less important as the dust self- 
heating is small compared to the dust heating from 
UV/optical photons. 

In practice, the placing of rays is controlled by the 
maximal numbers of rays per cell Ni and the "width" 
of the Tsou ~ 1 layers At; with |tsou - 1| < At/. The 
placing and iteration over wavelength is stopped when 
the change in the energy deposited in all cells drops 
below a chosen change limit and the mean field dis- 
tribution has been determined. If the energy does not 
converge after placing Ni rays in all cells within the 
Tsou ~ 1 layers, a second run with all rays is performed 
but using the updated energy information from the first 
run. If still no convergence is obtained in this run, Ni 
is increased to improve the resolution. The number of 
rays for typical 3D dust RT applications can easily ex- 
ceed 10^ 

The final step is to calculate the radiation that is re- 
ceived by the observer from each cell. This calculation 
makes use the pre-calculated Tobs- Moreover, radia- 
tion directly received from the inner and outer discrete 
sources as well as from the background radiation that 
is inside the field of view is calculated by placing rays 
from the sources and the background to the observer 

To illustrate how the rays are placed in various sit- 
uations, we present a few simple examples. The first 
is a molecular cloud core with a gas mass of 1 Mq 
that is illuminated by a strong MIR radiation field that 
dominates its thermal budget. The core has no in- 
ternal source of radiation except the thermally emit- 
ting dust. Furthermore, the self-heating of the dust 
by neighboring dust within the core can be neglected. 
The pre-calculations will not yield any T-surfaces since 
the optical depth is too low in the MIR. The main 
calculation will therefore be to propagate the external 
anisotropic radiation field through the core with little 
need for refinement by placing additional rays. Cor- 
respondingly, the rays can be placed dense enough in 



direction space to resolve all features of the extended 
illumination source (e.g., a nearby photodissociation 
region, see Steinacker et al. 2005). 

The second example is a binary star surrounded by 
a circumstellar disk. The pre-calculation will identify 
the inner disk and the disk surface as the t » 1-zone 
for the NIR wavelengths. The ray pattern will there- 
fore be dense at the inner rim, and in the zone above 
and below the disk. At MIR or longer wavelengths, the 
ray pattern will be less dense in the atmosphere since 
the T X 1 -range will move deeper into the disk as the 
thermally emitting inner dust is important. This appli- 
cation has an additional complication in that it requires 
the scattered light to be calculated (see below). 

More complex structures being illuminated by a 
central source are a problem for all RT solvers: a 
proper resolution requires many spatial cells and there- 
fore many rays from the star to the cell and from the 
cell into the surrounding regions. For RayT, if a fixed 
direction grid is used, more and more cells are over- 
looked when the diverging rays reach the outer parts. 
For such applications, beam-splitting can be used to 
split a ray in order to sample several neighboring cells 
(Abel & Wandelt 2002). 

Including scattered radiation 

While the thermal source contribution can be calcu- 
lated using the mean intensity J(x, A) per cell, the scat- 
tered light intensity depends on the direction of the in- 
coming radiation, the optical depth for scattering, and 
the phase function of the grains. The RayT scheme 
that includes thermal re-radiation and scattering must 
therefore store J and the intensity and direction of the 
incoming radiation for each cell and wavelength. 

Every RayT step is expanded to include "directly 
scattered light" by calculating the amount of radiation 
that is scattered towards the observer using the pre- 
calculated Tobs- In many applications, this singly scat- 
tered radiation is a good fraction of the total scattered 
light in the computational domain. Including multiple 
scattering in RayT methods significantly increases the 
computational requirements. Applications using 3D 
dust RT based on RayT with scattering sources have 
been presented for massive disk candidates (Steinacker 
et al. 2006) and grain growth in molecular cloud cores 
(Steinacker et al. 2010). Further codes using ray- 
tracing solvers and also alternative techniques to deal 
with high-optical depth are described in Pinte et al. 
(2009) and the references within. Fortunately, there 
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are specific astxophysical regimes where multiple scat- 
tering is not crucial or where reasonable approxima- 
tions can be made to limit the computational require- 
ments. 

The standard ISM dust grain distribution can be ap- 
proximated as a mean grain with a 0. 1 fim grain size. 
Such grains have very low scattering efficiency in the 
MIR and, thus, scattering is often neglected if calcu- 
lations are done at MIR or longer wavelengths. But 
caution should be exercised as for coagulated grains in 
dense molecular clouds scattering can be be important 
even at MIR wavelengths (Steinacker et al. 2010). In 
addition, dust grains have a strongly forward peaked 
scattering phase function and a portion of the scat- 
tered light can be transported along the already cal- 
culated rays (see, e.g., Steinacker et al. 2003). For the 
wavelengths where dust shows a more isotropic scat- 
tering phase function (e.g., MIR), additional rays can 
be added to carry the scattered radiation. An exam- 
ple of where adding scattering is fairly straightforward 
is in modeling star formation regions. Models of these 
regions will have a high number of rays in cells that are 
identified using the t ^ 1 -search algorithm, allowing 
scattering to be calculated along existing rays. Finally, 
the high albedo of dust grains and the distribution of 
the scattered intensity over the unit sphere naturally re- 
duce the intensity of the scattered radiation with sub- 
sequent scatterings allowing an intensity threshold to 
be used to limit the number of scatterings needing to 
be calculated. 

In astrophysical objects where the optical depth is 
very small, or very large, or where the t ^ 1 layer is 
not resolved, calculating multiple scattering can be ig- 
nored. If T <K 1 for an object (e.g., the diffuse ISM or 
the Zodiacal light), then single scattering completely 
dominates the scattered intensity and multiple scatter- 
ing can be ignored. In models of r » 1, the main mod- 
eling goal is usually to calculate the thermal dust emis- 
sion. Since scattering is generally important at optical 
and shorter wavelengths, most of energy at these wave- 
lengths is absorbed in the t x I regions of the object. 
Thus, not calculating the scattered light properly will 
not have a large impact on the dust emission results. 
For very embedded sources like massive young stellar 
objects, the stellar energy is completely converted to 
radiation at MIR or longer wavelengths before it leaves 
the core and therefore scattering does not influence the 
outer radiation field. Finally, if a model does not re- 
solve optically thick regions (e.g., in some extragalac- 
tic applications) and thus not not having t s; 1 regions. 



then ah the rays can be used to determine the global 
field and transport the single scattering. 

4.3. RayT Error analysis 

The precision of a solution of the RT equation along 
a ray through a spatial grid on which the opacity and 
source terms are described analytically is limited by 
machine precision. 

Solving the RT equation along a ray through a spa- 
tial grid on which the opacity and source terms are an- 
alytically described can be done with only machine 
precision errors. For example, the previously men- 
tioned Runge-Kutta solver with adaptive step size con- 
trol provides a good compromise between accuracy 
and computational cost while providing explicit error 
control. 

The main source of error for solutions on a sin- 
gle ray is the interpolation error of the density and 
source function from the underlying grid. It can accu- 
mulate if the density grid shows strong (and partially 
unresolved) gradients. The pre-calculation of t a; 1- 
regions and shielded areas helps to characterize how 
well the chosen grid describes the underlying physical 
problem, and refining the grid with this information 
can reduce the interpolation errors in the intensity dis- 
tribution. Unlike solvers calculating the intensity on a 
grid (e.g. using finite differencing), RayT solvers cre- 
ate no diffusive errors (beam-smearing). 

There is no general formula for the global error in 
the achieved intensity distribution, because deviations 
caused by coarse resolution due to the placement of 
rays are hard to quantify. A good test for the global ac- 
curacy of the thermal conversion of radiation is to cal- 
culate the energy output of the source term integrated 
over the domain, and to compare it to the energy in 
the radiation leaving the domain. Another test for the 
overall resolution of the important regions in the do- 
main is based on the number of crossing rays per cell. 
All the cells that contribute significantly to the radia- 
tion field overall should crossed by many rays. Having 
important cells with a small number of rays is an indi- 
cation of underresolving the grid. 

A practical test to understand and measure the 
global error is to increase the spatial or wavelength res- 
olution by inserting more rays or adding wavelength 
grid points, and testing the stability of the solution. 
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5. THE MONTE CARLO SOLUTION METHOD 

The Monte Carlo (MC) method is a general compu- 
tational technique that is widely used in many different 
areas, including numerical mathematics, physical sci- 
ences, finance, and medicine. It is a joint name for a 
variety of stochastic or probabilistic techniques, which 
all have in common that they solve equations by sam- 
pling random numbers. MC methods are particularly 
interesting for transport systems, which was the moti- 
vation of their first application in the 1940s. For a gen- 
eral overview of MC as a tool for transport problems, 
see e.g. Dupree & Fraley (2002), Kalos & Whitlock 
(2009) or Whitney (2011). Its application to dust RT 
problems in an astrophysical context has a history of 
more than 40 years (e.g Mattila 1970; Roark, Roark 
& Collins 1974; Witt & Stephens 1974; Witt 1977). 
In the last four decades, it has become a mainstream 
method for 3D dust RT calculations. 

The basis of MC RT is to treat the radiation field as 
the flow of a large but finite number of photon pack- 
ages (often called photons). Each individual photon is 
followed along its journey through the dusty medium. 
At every stage in its journey, the characteristics that de- 
termine the path of each photon (such as its birth loca- 
tion, initial propagation direction, or the distance along 
the path until the next interaction with a dust grain) are 
determined in a probabilistic way by generating ran- 
dom numbers from an appropriate probability density 
functions (PDF). At the end of the simulation, the ra- 
diation field is recovered from a statistical analysis of 
the photon paths. Hence, the MC technique simulates 
the RT instead of explicitly solving the RT equation. 

Central to all MC techniques is the process of gen- 
erating random numbers from a given PDF p(x)dx. 
Thus an algorithm is needed that returns a set of num- 
bers X such that the probability that X lies in the in- 
finitesimal interval between x and x + dx is equal to 
p{x) dx. The starting point for such algorithms, which 
are key to the MC process, is a (pseudo-)random num- 
ber generator. This is a code that generates uniform 
deviates (random numbers with an equal probability to 
be chosen in the unit interval between and 1). In 
order to generate random numbers from another, arbi- 
trary PDF, one almost always applies appropriate op- 
erations on one or more uniform deviates. The most 
popular methods are the so-called transformation and 
rejection methods, details of which can be found in 
Devroye (1986) or Press et al. (2002, Ch. 7). 

In the following subsections, we will gradually in- 



troduce the ingredients and techniques that are com- 
bined to to develop a 3D MC dust RT code. We start 
by describing the simplest MC RT technique, as it is 
the basis for all modern MC RT, in §5.1. These sim- 
ple techniques are sufficient for MC calculations in 
geometries with a large degree of symmetry, such as 
ID spherical or slab geometries. For 3D applications, 
however, they would result in codes that would be very 
inefficient. Fortunately, there are various weighting 
schemes that improve the performance of this tech- 
nique making modern MC RT quite efficient. Some 
of these weighting schemes have a strong heritage in 
RayT methods, making most modem MC RT codes 
hybrids between the classical MC and RayT tech- 
niques. A number of these weighting schemes were 
developed for 2D geometries, especially cylindrical 
geometries. The use of weighting schemes is critical 
for 3D MC RT; in §5.2 we discuss several of these 
techniques. 

5.1. Simple MC RT 

The simplest MC calculation consists of consider- 
ing the RT problem at a single wavelength A. We con- 
sider a source of photons, characterized by the source 
term jt(x, n. A), and a distribution of dust, character- 
ized by the dust density p{x). Throughout the cal- 
culation, the state of each photon is tabulated by its 
energy, position, direction of travel, and polarization 
state. For 3D RT, the Cartesian coordinate system is 
usually used resulting in the photon's position being 
stored as jc = {x,y,z) and the direction using direction 
cosines as « = (rix, riy, n^). The polarization state is 
stored using the Stokes vector, S - {I, Q, U, V). 

Step 1: birth 

The first step in a photon's life cycle is its birth, i.e. 
its injection into the computational domain. If is the 
number of photons in the model run and Liot(A) is the 
total luminosity of the source, the luminosity carried 
by each photon is L = Ltot(A)/N. The initial position x 
and propagation direction n are to be chosen randomly 
according to the source term jt(x,n,A), which means 
that they need to be sampled from the PDF 



p(x, n) dx dn 



jt(x,n,A)dxdn jjx,n. A) dx dn 



jj jt{x, n, A) dx dn 



itot(^) 



(23) 

In many cases, e.g. for emission by stars or thermal 
emission by dust grains, the emission is isotropic, 
which implies that the initial propagation direction can 
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be chosen randomly from the unit sphere, 

dn sinOdOdd) 
/7(«)dn= — = -. (24) 

4-7T 4-7T 

Generating random values for 6 and (p from this dis- 
tribution can easily be done using the transforma- 
tion method. In other cases, e.g. when we deal with 
external emission illuminating an interstellar cloud 
or anisotropic emission from an AGN accretion disc, 
different PDFs for the initial propagation direction 
need to be considered (e.g. Niccolini, Woitke & Lopez 
2003; Stalevski et al. 2012). For most cases, the pho- 
ton is assumed to be emitted unpolarized, i.e. S = 
(1,0,0,0). 

Step 2: determination of the interaction point 

Once the photon is launched into the dusty medium, 
the next step consists of randomly determining whether 
it will interact with a dust grain, and if so, where this 
interaction will take place. The PDF that describes the 
free pathlength before an interaction is most conve- 
niently described in optical depth space, where it has 
an exponential distribution, p{t) dr = e^'^dr. The opti- 
cal depth T to which a particular photon travels along 
its path before it interacts with a dust grain, is drawn 
from from this exponential distribution. This is easily 
done using the transformation method: we simply pick 
a uniform deviate ^ and solve the equation 

^= r e-^'dr' (25) 
Jo 

for T. Integrating and solving gives r = - In ^ re- 
membering that the distributions of ^ and 1 - ^ are 
equivalent. If r is greater than the optical depth Tpath 
to the surface of the system in the direction the pho- 
ton is traveling, the photon escapes the system, and the 
hfe cycle of this particular photon is over Otherwise, 
the photon will interact with the dust medium at the 
location along the path corresponding to the traveled 
optical depth t. 

The next step is converting the traveled optical 
depth T to a physical path length s, such that we can 
move the photon to the interaction site. This means 
that we have to integrate along the path and solve the 
integral equation 

I Ken(s',yl)p(s')ds' ^ T (26) 
Jo 

for the path length s. Comparison between this equa- 
tion and Eq. (22) highlights the intimate link between 



the MC solution technique and the RayT technique: a 
large fraction of the calculations in MC simulations are 
pure ray-tracing operations. 

In practice, MC codes virtually always use a grid 
structure of dust cells on which the dust density and 
the optical properties are discretized. The integral in 
equation (26) then reduces to a sum over the consecu- 
tive grid cells along the path, and the integral equation 
comes down to summing the optical depth along the 
path until t is reached. This calculation is often one 
of the more computationally intensive portions of MC 
RT and is a strong driver for choosing a grid optimized 
for the particular astrophysical object of interest. For a 
Cartesian model grid, the distance traveled in each cell 
is easy to calculate, as is the next grid cell along the 
path. For hierarchical grids, more complex neighbor 
search algorithms may be required. 

Step 3: absorption and scattering 

Once the path length s has been calculated, the pho- 
ton moves from its old location x to its updated loca- 
tion, i.e. the interaction site x + s n. At this location, 
the photon can either be absorbed or scattered; the ap- 
propriate PDF is hence not a continuous, but a discrete 
function with only two possible values. The probabil- 
ity that the interaction is a scattering event is equal to 
the dust albedo a - /Csca/'<'ext- Using a uniform deviate 
^, the nature of the interaction is easily determined: if 
^ < fl we have a scattering event, if ^ > a an absorption 
event. 

In the case of an absorption event, this is the end of 
the photon's life cycle. If dust emission is to be cal- 
culated in the simulation, the absorbed photon lumi- 
nosity is stored in the interaction cell. This absorbed 
luminosity will be used at a later stage to calculate the 
dust emission spectrum, which can then be used as the 
source function for another MC cycle. 

If the interaction is a scattering event, the next step 
consists of determining the new propagation direction. 
In the case of isotropic scattering, this just comes down 
to generating a new random point from the unit sphere. 
In the case of anisotropic scattering, the new propaga- 
tion direction n is chosen according to the PDF 

Mn,n',x,A)dn 

p{n) dn = , (27) 

4-n 

where <S>{n,n',x,A) is the scattering phase function 
and «' is the original propagation direction before the 
scattering event. For spherical dust grains, the scat- 
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tering phase function is dependent only on the scat- 
tering angle between n and n' . For the HG phase 
function, the most popular approximation to the real 
phase function, the generation of a random scattering 
angle and hence the calculation of the new propaga- 
tion direction, can be done analytically (Witt 1977). 
For polarized RT, the scattering process is more com- 
plex, as the scattering phase function is dependent on 
the polarization of the photon and, in addition, scat- 
tering changes the polarization state of the photon. A 
detailed description of scattering events in the case 
of polarized radiation can be found in e.g. Fischer, 
Henning & Yorke (1994), Code & Whitney (1995), or 
Goosmann & Gaskell (2007). 

With its new propagation direction determined, the 
photon can continue its journey through the dusty 
medium. This means that the second and third step 
are repeated until the photon either escapes from the 
system, or is absorbed by a dust grain. ^ This entire 
process is repeated for all the photons until the last 
photon has left the dusty medium. The RT at differ- 
ent wavelengths can be independently calculated given 
there is no explicit wavelength coupling in the RT. 
There is an imphcit coupling when dust emission is 
included in the modeling. In this case, the absorbed 
luminosity at every wavelength is stored in every dust 
cell in the spatial grid. After finishing the simulation 
at all the wavelengths, the absorbed luminosity is used 
to compute the mean intensity J{x, A) of the radiation 
field and subsequently the dust source term jiix, A). It 
is common to include the dust emission as a second 
source of photons by first computing the RT from the 
primary sources, computing the RT for the dust emis- 
sion, recomputing the dust emission with the new ra- 
diation field, and iterating until a set convergence level 
is reached. 

In simple MC RT, the output desired from the calcu- 
lation is usually the view of the system from a partic- 
ular observer location. Images of the system are con- 
structed by projecting the position of each photon that 
escapes in the direction of the observer within some 
angular tolerance onto the plane of the sky. This plane 



*It is possible to have cases wliere the optical depth is veiy large and 
the number of scatterings is coiTespondingly large. It is common to 
impose a limit on the number of scatterings to calculate as a result. 
This has the consequence that high optical depths are not sampled 
well with potential systematic errors in the calculation. For most 
cases, the systematic errors from imposing a maximum number of 
scatterings are negligible. 



is divided into pixels and the images are built from 
those photons. A large number of the photons will es- 
cape the system in directions other than the observer 
and are not counted, unless symmetries (e.g., circular 
or cylindrical) can be exploited. Special care must be 
applied when constructing the images which give the 
polarization state of the scattered flux, specifically the 
Q and U images. As part of the construction of these 
images, the polarization vector of the photon must be 
rotated so that it is referenced correctly in the image 
plane. 

5.2. Weighted MC RT 

Simple MC RT is the easiest to understand, but it 
would be extremely inefficient for full 3D calculations. 
For ID and 2D MC RT calculations, there are sym- 
metries, such as spherical and cylindrical, that can be 
exploited and thus the number of photons needed to 
achieve accurate results is relatively small. In the case 
of 3D dust RT, there are no symmetries by definition; 
this has been one of the motivations for a number of ac- 
celeration methods. Most of these acceleration meth- 
ods are well established and validated, while others 
still have a more experimental character. Almost all 
of these methods were first developed for 2D RT cal- 
culations. The benefits of weighted MC techniques as 
compared to simple MC RT methods are illustrated in 
Fig. 4. 

The basis of all acceleration methods is to assign 
a weight W to each photon and modify this weight. 
The weight of each photon is equivalent to the frac- 
tion of the luminosity of the emission source carried by 
that photon i.e. the number of photons in each photon 
package. Several acceleration techniques use the idea 
of biasing, i.e. generating random numbers from a PDF 
q{x) Ax rather than from the appropriate PDF p(x) dx. 
This biased behavior is corrected for by assigning the 
weight W - p{x)/q(x) to the photon. The biasing tech- 
nique is used in many MC applications, and can be a 
very effective way of reducing the variance (Dupree & 
Fraley 2002). It is, for example, common practice in 
MC numerical integration where it is known as impor- 
tance sampling. 

Biased emission 

A direct application of the biasing technique is the 
so-called biased emission. The initial propagation di- 
rection of the photons launched into the dust medium 
is usually determined from the angular part of the 
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source term j{x, n. A). There are cases, however, where 
an increased level of emission in particular directions 
is desired, for example to increase the signal-to-noise 
in particularly interesting directions, such as the polar 
regions of a star with an accretion disk (Yusef-Zadeh, 
Morris & White 1984). In this case, the emission of 
photons is biased towards the directions of interest, 
and the initial weight of the photon is determined us- 
ing the standard biasing weight formula. The same 
technique can also be applied to the spatial part of the 
source term to increase the number of photons emitted 
from regions with a low emission rate (Juvela 2005). 

This technique of biased emission has the potential 
of strongly increasing the efficiency of a MC simula- 
tion. On the other hand, it has the drawback that it is 
very model-dependent, and therefore requires signifi- 
cant manual interaction. It is, in a sense, comparable 
to the placement of the rays in the RayT method. 

Absorption-scattering split 

This acceleration method allows for a photon to 
contribute to both absorption and scattering at each in- 
teraction site. Instead of randomly choosing the nature 
of the interaction, the photon is split in two parts: one 
that is absorbed and one that scatters. The fraction ab- 
sorbed is equal to (1 - a) times the current weight of 
the photon. For the scattered part that continues its life 
cycle through the dust, the weight is multiplied by a 
factor a. 

Forced scattering 

Instead of having each photon either scatter or es- 
cape the system, the scattering can be forced to occur 
every time (Cashwell & Everett 1959). In the simple 
MC routine, the randomly generated optical depth is 
compared to the total optical Tp^th along the photon's 
path. This approach has a problem for regions with 
low optical depth: in those regions, many photons just 
leave the system without interacting with the dust, re- 
sulting in a low efficiency of dust scattering and heat- 
ing. A way to avoid this low efficiency is the tech- 
nique of forced scattering, which limits the values of 
the randomly chosen optical depths to the range be- 
tween and Tpath- This can be achieved by biasing 
the PDF from which the optical depth is generated. In- 
stead of sampling from the actual exponential PDF, we 
consider an exponential distribution cut off at t = Tpath- 



Properly normalized, this PDF reads 

( e-^ dr 

g(T)dT = i 1 - e-^p'»i. ^<''"path, 1^28) 

[ T> Tpath. 

Generating a random t from this distribution is 
straightforward, and guaranteed to produce an inter- 
action before the photon has left the system. The com- 
pensation to be applied to the weight of the photon is 
a factor 

Wft = ^ = 1 - e-^P"'" . (29) 
P(t) 

One issue with forcing every scattering, when com- 
bined with the absorption-scattering split, is that there 
is no longer a natural stopping criterion for the scat- 
tering calculation. In the original MC cycle, photons 
end their journey when they are either absorbed by the 
dust or leave the dusty system. The common solution 
is to set a minimum weight for a photon, below which 
the photon's life cycle is terminated. The value of this 
termination weight is usually set to be very low, after 
the equivalent of many scatterings. An alternative so- 
lution is to only force the first scattering and revert to 
the standard scattering calculation for subsequent scat- 
terings. 

Peel-off technique 

For 3D cases, it is usually desired to calculate the 
appearance of a system for an observer at a partic- 
ular orientation. Simple Monte Carlo RT is particu- 
larly inefficient in building up such an image as only 
the photons that are emitted from the system in the 
direction of the observer contribute to the output ap- 
pearance. In addition, some blurring of the image is 
inherent as photons that are within some tolerance of 
the desired direction are used. This inefficiency can be 
eliminated by requiring that all photons directly con- 
tribute to the output images by calculating the portion 
of the photon that is emitted from sources and scat- 
tered at every interaction point in the observer's direc- 
tion (Yusef-Zadeh, Morris & White 1984). The weight 
factor of a photon in the direction of the observer is 

H'po =p(«obs)e-"* (30) 

where Tobs is the optical depth from the position of 
the emission or scattering event, and piiiob^) is the 
probability that the photon would be directed towards 
the observer. For example, for isotropic emission, 
pitobs) - 1; and after a scattering event we have 
piiiabs) = ^in,nobs,x,A). 
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Peel-off is a pure application of RayT, again high- 
lighting the connection between both approaches. 
Peel-off is probably the most important acceleration 
technique for 3D MC RT simulation. It has a sig- 
nificant computational cost, however, as it requires a 
calculation of the optical depth from the emission or 
scattering location towards the observer after every 
emission or scattering event. One way to alleviate this 
computational burden is by precalculating the optical 
depth towards the observer for each cell. Another pos- 
sibility is to store the information per cell and create 
the images at the end of the simulation (e.g. Dulle- 
mond & Turolla 2000; Pinte et al. 2006; Min et al. 
2009). Using a precomputed r^hs or using stored infor- 
mation to compute the image does come with a price, 
loss of subgrid resolution in the model images. The 
benefit of this approximation in decreased computa- 
tion time has to be weighted against the loss of subgrid 
resolution for the particular modeled astrophysical ob- 
ject. 

Continuous absorption 

The accuracy of the dust re-emission of the ab- 
sorbed energy can be enhanced by absorbing not just 
at the interaction site, but along the path the photon 
travels. Depending on the implementation, this can be 
done up to the location of the scattering event (Lucy 
1999), or along the entire path (Niccolini, Woitke & 
Lopez 2003; Baes et al. 2011). In the latter scenario, 
the photon is effectively split in -H 2 different parts: 
one part Wesc that leaves the system (and is hence not 
accounted for anymore), one part Wsca that is scattered 
at the interaction location (determined stochastically), 
and parts Wabsj that represent the fraction absorbed 
in each of the cells along the photon's path. These 
different fractions are 

Wesc = e"'^"'" , (31) 
Wsca = fl(l-e-^^'"^), (32) 
W,i,,j^(l-a)ie-''-' -e-''), (33) 

where tj is the optical depth measured from the pho- 
ton's location to the surface of the j'th cell along the 
path. An alternative interpretation of this continuous 
absorption approach is that it estimates the mean inten- 
sity of the radiation field not by means of the absorbed 
luminosity, but by means of counting the luminosity 
that passes through each cell. The strength of the tech- 
nique lies in the fact that all photons contribute to the 
calculation of the absorption rate of each cell they pass 



through, and not only of those cells with which they in- 
teract. This is particularly useful for the optically thin 
regime, which has very few absorptions in the simple 
MC approach. 

Instantaneous dust emission 

The traditional method in computing the self- 
consistent dust emission in a MC RT simulation con- 
sists of running an independent MC simulation at ev- 
ery individual wavelength and storing the absorbed 
luminosity over the computational grid. In a sec- 
ond stage, the dust emission spectrum is calculated 
in each dust cell and used as secondary source term. 
This approach inevitably leads to iteration, as the dust 
emission itself affects the radiation field. It is pos- 
sible to compute the output dust emission spectrum 
from a multi wavelength model without iterating. The 
instantaneous dust emission technique or frequency 
distribution adjustment technique emits dust thermal 
photons immediately after each absorption event in a 
specific grid cell, with the wavelength of the emitted 
thermal dust photon carefully chosen to retain ther- 
mal equilibrium (Bjorkman & Wood 2001, Baes et al. 
2005). 

Besides eliminating the need for iteration in the 
computation of the dust emission spectrum, the instan- 
taneous dust emission technique also has the advan- 
tage that photon packages are emitted at the exact po- 
sition where the absorption event took place, so sub- 
grid resolution is achieved. One disadvantage is that 
the absorption/re-emission event takes place in a sin- 
gle cell, which results in a poor convergence rate of 
the radiation field in cells with a low dust absorption 
rate. This makes the classical iterative technique with 
continuous absorption more efficient than the instanta- 
neous dust emission technique, at least when applied in 
its original form, for 3D simulations (e.g. Chakrabarti 
& Whitney 2009). This problem can be alleviated 
by applying a combination of the instantaneous dust 
emission and the continuous absorption techniques. In 
this hybrid method, the photon packages are followed 
through the domain using the instantaneous dust emis- 
sion technique, but the final dust emission spectrum 
of the cells, used to create images and SEDs, is cal- 
culated based on the continuous absorption approach 
(Pinte et al. 2006). A second problem for the instan- 
taneous dust emission technique is that it was origi- 
nally designed to work with equilibrium dust emission. 
Several ways have been explored to adapt the instan- 
taneous dust emission technique for transiently heated 
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grains (Kriigel 2007, Wood et al. 2008, Heymann & 
Siebenmorgen 2012). 

High optical depths 

In regions of high optical depth, the simple MC 
routine becomes very inefficient, as photons can be 
trapped in a virtually never-ending loop of scatter- 
ing events. This problem is largely solved when the 
absorption-scattering split is applied, as the weight of 
the photon then decreases with every scattering event 
terminating when a very low weight is reached. How- 
ever, in regions with extreme optical depths (such as in 
the midplane of circumstellar discs around protostars) 
or at wavelengths where scattering largely dominates 
absorption (such as the far-UV), this can still imply 
a significant computation burden. A solution to this 
problem is to mirror or reflect photons from high opti- 
cal depth regions and use the diffusion approximation 
to find the RT solution in these regions. The diffusion 
approximation allows for multiple interaction steps to 
be calculated in a single computation. An elegant solu- 
tion that is well adapted to the MC method is to solve 
for the RT in optically thick cells using a modified ran- 
dom walk technique that also uses the diffusion ap- 
proximation (Min et al. 2009; Robitaille 2010). 

Polychromatism 

For multi-wavelength RT, it is possible to signifi- 
cantly speed up the calculation by considering photon 
packages that consist of photons of all wavelengths. 
The advantage of this technique that a MC run is si- 
multaneously solved at all wavelengths, instead of a 
run for each wavelength. The difficulty in this ap- 
proach is that many of the PDFs that describe the 
Ufe cycle of a photon are wavelength dependent, such 
as the path length distribution or the scattering phase 
function. One solution is to consider partly polychro- 
matic photon packages, which shift to monochromati- 
cism as soon as wavelength-dependent PDFs are in- 
volved (Baes, Dejonghe & Davies 2005). A more ad- 
vanced option is perform the calculations at one ref- 
erence wavelength /i,ef and use the biasing technique 
to adjust for the wavelength-dependent PDFs (Juvela 
2005; Jonsson 2006). The efficiency gains of full poly- 
chromatic RT are large given that every photon calcu- 
lated would contribute to the output images and radia- 
tion field density at all wavelengths instead of a single 
wavelength. But there is a known significant compli- 
cation - the biasing factors can be very large and, as 



a result, can dominate the results at a particular wave- 
length. This has a systematic effect on the results that 
is hard to control and therefore, this method should be 
considered experimental. 

5.3. Uncertainties for Monte Carlo 

In the simple (unweighted) MC RT, the noise in the 
output quantities (i.e., scattered intensity, polarization, 
etc.) scales as N^^^^, where N is the number of pho- 
tons. In the case of weighted Monte Carlo, the uncer- 
tainties in the output quantities does not scale directly 
withA^-i/2_ 

The uncertainties can be calculated by using the dis- 
persion in the average properties of the photons used to 
determine an integrated quantity (Gordon et al. 2001). 
If the integrated quantity is X, then 

N 

X = Y^xj ^ NJ, (34) 

where xj is the contribution of the /th photon to X, N 
is the total number of photons in the model run, and 
X is the average contribution each photon makes to X. 
The uncertainty in X is then crx - X {cTxIx), where cr ^ 
is the standard deviation of x, calculated using 

The equations for the uncertainties in output quantities 
given above can be used to explicitly enforce a par- 
ticular uncertainty level in the final results of a model 
run. The uncertainties in the output quantities of in- 
terest can be checked during the model run and the 
number of photons adjusted dynamically to achieve 
the desired accuracy. For example, when calculating 
a multi-wavelength model more photons can be put at 
the wavelengths with higher optical depths compared 
to wavelengths with lower optical depths. As the num- 
ber of photons is determined from the statistics of the 
particular run itself, it automatically takes into account 
the full RT solution, including the locations of the pho- 
ton sources, dust, etc. The output quantities to be used 
for uncertainty control can range from the global flux, 
to pixels in resolved images, to the radiation field den- 
sity in each cell. Additionally, it is possible to improve 
the convergence to the needed accuracy by identify- 
ing areas of the model with high noise (e.g., partic- 
ular cells with few photon interactions) and dynami- 
cally adding photons to those areas (e.g., though bi- 
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ased emission). The use of dynamically determine un- 
certainties to improve the convergence of a model has 
been done in limited cases and is clearly an area for 
future improvement. 

Uncertainties associated with the model setup itself 
are more difficult to quantify. Evaluating the system- 
atic uncertainties due to specific model choices is usu- 
ally done by trying other parametrizations or increas- 
ing the grid resolution and evaluating how the out- 
put quantities change. Example model setup choices 
that are prone to systematic uncertainties are the dust 
density grid, specific parametrizations of the photon 
sources, and wavelength grid. These kinds of uncer- 
tainties are clearly the most diflicult to diagnose and 
generally rely on the expertise of the coder and user. 

6. CHALLENGES IN MODELING OBSERVA- 
TIONS 

The uses of 3D dust RT codes are many, from 
understanding the impact of locally clumpy dust on 
dust RT (Witt & Gordon 1996, 2000; Bianchi 2008) 
to modeling observed images of objects to derive the 
source and dust distributions (De Looze et al. 2012b; 
Steinacker et al. 2005; Gordon et al. 1994) to predict- 
ing the appearance of an object that has been modeled 
with a HD code (Steinacker et al. 2004; Bethell et al. 
2004; Jonsson, Groves & Cox 2010; Robitaille 2011) 
to investigating the ensemble behavior of objects to de- 
rive average source and dust properties (Law, Gordon 
& Misselt 201 1) to test observability of objects for in- 
strument and observation planning (Wolf 2008; Gon- 
zalez et al. 2012). To illustrate the challenges in mod- 
eling observations, we focus on the modeling of ob- 
jects to derive their physical properties from observa- 
tional data. The common approach is to perform some 
steps of the following scheme: 

1 . Choose a parametrized model for the dust den- 
sity, radiation sources, and dust optical proper- 
ties. 

2. Discretize the problem with the choice of grids 
and use the RT code to derive a SED and/or im- 
ages for this object model. 

3. Compare simulated and observed SED or im- 
ages. 

4. Evaluate the differences in the SED or images 
and change the model parameters and/or the 



model assumptions to minimized these dififer- 
ences. Repeat step 2. 

5. Find the parameter sets that come closest to de- 
scribing the observed data. 

6. Evaluate the observational and theoretical errors 
that are important for the comparison. 

Although the scheme seems straightforward, the full 
set of steps has rarely been carried out for existing, 
published 3D dust RT models as each of these steps 
provides challenges for 3D dust RT. 

Model choice 

The aim is to choose an appropriate model that al- 
lows meaningful physical information to be derived 
from the given observational data. This is a general 
topic, but 3D RT modeling has special features which 
are important to consider The data are usually SEDs 
and/or images and we can investigate the number of 
independent information bits a priori, e.g. by counting 
the number of wavelengths for which the flux density 
values have been measured. This is then compared to 
the number of model parameters. 

Most 3D models include several 10 to hundreds of 
free parameters depending on the complexity of the 
spatial model and the assumed dust properties. Even 
for attempting only to reproduce the global SED of 
an object, it is difficult to stay below 10 parameters 
since modeling a dusty 3D structure requires a few 
length scales and/or power laws to describe density, 
dust properties, and viewing angles. Fig. 5 illustrates 
the large number of free parameters needed to repro- 
duce the 3D dust distribution of a molecular cloud. 
Fig. 6 shows four examples of state-of-the-art RT fits 
to dusty galaxy images. Such RT model fits are de- 
signed to determine the parameters that describe the 
intrinsic 3D distribution of stars and dust both large- 
scale geometric parameters (e.g., stellar and dust scale 
lengths and heights) and measures of the small scale 
inhomogenity or dumpiness of the dust. 

Starting with more free model parameters than data 
points is often considered to produce meaningless pa- 
rameter values. However, a careful exploration of the 
parameter space can detect ambiguities and reveal re- 
dundant parameters. As long as the parameter space 
exploration can be afforded computationally, starting 
with a detailed model and analyzing the parameter am- 
biguities is usually the most accurate way to model an 
object. 
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Gridding 

As discussed in §3, there are a variety of options 
for discretizing space, direction, wavelength, and the 
dust model. If the grid fails to resolve structures that 
are important at a particular wavelength, the corre- 
sponding model image will be correspondingly inac- 
curate. A good example is the surface layer of an 
accretion disk illuminated by a central star. If the 
layer is not well-resolved, the description of the ra- 
diation field that is scattered in the layer will be poor, 
and the scattered light images inaccurate. Therefore, 
the discretization step often includes running the RT 
code before the start of the modeling to verify that ex- 
pected features are present in the image or SED, and 
that changes/refinement do not influence the overall 
appearance of the object. 

Comparison of model and data 

Once the model images or SED have been calcu- 
lated, the results should be convolved with the beam 
of the observed instruments/telescopes and the sam- 
pling should be made equal (pixel size for images). It 
is important to apply detection limits, especially when 
studying faint structures. For interferometric obser- 
vations, the incomplete coverage of all spatial scales 
implies that the comparison should be performed in 
the (u,v) plane rather than in spatial plane to achieve 
the highest fidelity. Calibration uncertainties and cor- 
related noise properties in the observations should be 
understood and included in the model versus obser- 
vations comparison. In addition, foreground emission 
(e.g. caused by the zodiacal light) should be carefully 
removed or modeled, especially for on-off" chopping 
observing modes that remove this emission as part of 
the data reduction process. Generally, modeling can 
benefit from good communication with experienced 
observers. Another potential source of error in com- 
paring model results and data can be a displacement 
of structures due to observational uncertainties. Addi- 
tional translation parameters can correct this and pro- 
vide physical insight by deriving improved position- 
ing. 

Exploration of the parameter space 

The parameter space of 3D dust distributions is 
large: a uniform coverage of a 10-dimensional pa- 
rameter space with 5 grid points in each parameter, 
would imply about 10 million individual RT calcula- 
tions. Therefore, almost all 3D RT modeling of data 



has been done "by hand," that is starting from a point 
in the parameter space and then varying just one pa- 
rameter to explore the variations in the results. This 
drastically reduces the coverage in parameter space at 
the expense of not exploring correlations between pa- 
rameters. In structures with strong gradients, inhomo- 
geneously distributed radiation sources or varying op- 
tical depth, the resulting radiation field may strongly 
react to changes in the model parameter. In this case, 
low-coverage "by hand" explorations are likely quite 
unreliable. 

To evaluate the model and the variation in the model 
parameter, the difiference between the model data and 
the observed data is defined (often using a metric) 
and the fitting procedure is then to minimize this dif- 
ference by varying the parameters. There are standard 
optimization methods that have been applied in astro- 
physics such as gradient descent, Newton's method. 
Metropolis optimization, or genetic algorithms. The 
latter two are able to leave local minima with the goal 
of ending in the deepest minimum providing the best- 
fitting model parameters. 

There are advancements in the application of au- 
tomated fitting techniques that are already in use for 
2D RT calculations and are or should be used for 3D 
RT calculations as well. One such technique is the 
Metropolis algorithm that has been applied to 3D dust 
RT in the form of simulated annealing (Steinacker 
et al. 2005). Precomputing a large grid of models is 
a promising technique. Robitaille et al. (2007) used 
a parametrized circumstellar disk model to create just 
such a large grid of precalculated SEDs for a 2D con- 
figuration and then used the grid to fit the disk parame- 
ters and to characterize uncertainties in fit parameters. 
Other promising techniques are 2D RT fitting tech- 
niques based on the Levenberg-Marquardt algorithm 
(e.g. Xilouris et al. 1999), downhill-simplex method 
(Bianchi 2007), and genetic algorithms (De Geyter 
et al. 2012). Some of these techniques have already 
been applied to 3D structures, albeit in limited param- 
eter spaces (e.g. Witt & Gordon 2000; Schechtman- 
Rook, Bershady & Wood 2012). 

Quantifying the ambiguity of the derived parame- 
ters is also important. The 2D circumstellar disk mod- 
eling of SEDs in Robitaille et al. (2007) provides a 
template for determining the ambiguity when a com- 
plex model is applied to only a few SED points. A 
general strategy to assess ambiguity is to explore the 
variation of the fitting metric (e.g., x^) in the vicin- 
ity of the best fit parameters, as well as to characterize 
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the overall degree of variation by random or grid-based 
parameter space exploration. 

Error analysis 

Each of the RT solvers has its own sources of errors 
due to approximations, the underlying grid, undetected 
regions of the computational domains, etc. Tracing 
these errors, measuring them, and estimating their im- 
portance is non-trivial. RayT provides error control for 
the solution along one ray and MC can provide error 
control on the number of photons needed, but whether 
the rays are sent through the right parts of the domain 
is not well quantified. For modeling, a comparison of 
these errors with the solution and observational errors 
determines the sensitivity of the observed data to the 
specific science question of interest. 

Inverse RT 

The true challenge of RT modeling is to invert it: 
determine the 3D density structure of the dust and 
the sources, and the dust properties from a set of im- 
ages taken at different wavelengths. Direct inversion 
is computationally impossible in 3D with current or 
expected computing capabilities. In ID an analytic in- 
version can be performed under special assumptions 
(Steinacker, Michel & Bacmann 2002), and numeri- 
cal ID inverse RT modeling was used in Doty et al. 
(2010). A forward RT method using many iterations 
provides a good solution, but is limited to fairly sim- 
ple RT applications. A prominent difficulty is the loss 
of information due to the line-of-sight integration in- 
herent in the observations. Multiwavelength observa- 
tions can help disentangle the line-of-sight integration. 
For example, points in the center of a molecular core 
are better shielded and cold compared to points in the 
outer parts, so mm observations can be used to con- 
strain the core center and shorter wavelength observa- 
tions the outer parts. Such molecular cores are simple 
enough to enable inverse RT to be done. For example, 
Steinacker et al. (2005) used a 3D background radia- 
tion field illuminating a 3D model core and fitted the 
MIR and mm images of the core Rho Oph D, deriv- 
ing the density and temperature structure with assumed 
dust properties. 



7. CODES AND BENCHMARKS 

7.1. Available 3D codes 

Until the mid-1990s, few RT codes could handle the 
absorption, scattering, and thermal emission by inter- 
stellar dust in general 3D geometries. The available 
codes were either limited to ID or 2D geometries, or 
were not able to calculate full 3D problem (e.g., miss- 
ing either scattering or thermal emission). The spec- 
tacular increase in computational capabilities during 
the past two decades, as well as the development of 
more powerful techniques to solve the RT problem, 
has led to the creation of a number of 3D RT codes. 
There are now (to our knowledge) almost 30 codes op- 
erational that can handle the full dust RT problem, i.e. 
including absorption and multiple anisotropic scatter- 
ing, in a general 3D geometry. Table 1 gives a list 
of published 3D dust RT codes that have been or are 
being used for astronomical applications. The applica- 
tion fields of the different codes varies widely, ranging 
from prestellar cores and circumstellar discs to active 
galactic nuclei and galaxies. 

With many different codes and several techniques 
available, it might be difficult for a potential user or 
future developer to decide which one to choose. The 
choice for a given code or solution technique should 
primarily be driven by the specific nature of the prob- 
lem that is to solved. While most of the 3D codes in 
Table 1 are applicable for a wide range of RT prob- 
lems, virtually all of them been developed with a par- 
ticular application in mind, and hence have been opti- 
mized for that particular goal. 

The programming language can also influence the 
choice; most existing 3D RT codes have been coded 
in FORTRAN, C, or C++. There are no features in 
the 3D RT problem itself or the two solution methods 
presented in this review that make certain languages 
preferential to others. The main driver is speed: since 
the the 3D RT problem is computationally very chal- 
lenging, all codes should be developed in a language 
suitable for high-performance computing. Moreover, 
parallelization is becoming increasingly important as 
a means to increase computational speed and mem- 
ory availability. Several codes are designed to work on 
shared-memory or distributed-memory clusters and/or 
adopt graphical processing units (e.g. Jonsson 2006; 
Jonsson & Primack 2010; Baes et al. 2011; Robitaille 
2011; Heymann & Siebenmorgen 2012). 

The choice of a given code or method primarily de- 



24 



Table 1: List of published 3D dust RT codes that have been or are being used for astronomical applications. The codes 
are ordered alphabetically as it is difficult to establish true code development dates given that many codes are first used 
for specific investigations with the full code details (references in table) published at a later date. 



name 


type" 


reference 


main application 


SKIRT 


MC 


Baes+ 2003; 201 1 


galaxies, AGNs 


(no name) 


MC 


BetheII+ 2004; 2007 


SE clouds 


TRADING 


MC 


Bianchi+ 1996, Bianchi 2008 


galaxies 


RADISHE 


MC 


Chakrabarti+ 2007; 2009 


galaxies 


(no name) 


MC 


Doty+ 2005 


SE clouds 


RADMC-3D 


MC 


Dullemond (in prep.) 


SE disks 


MOCASSIN 


MC 


ErcoIano+ 2005 


photoionized regions 


(no name) 


MC 


Eischer+ 1994 


SE disks 


(no name) 


MC 


Gon{aIves+ 2004 


SE clouds 




MC 


Goosmann+ 2007 


AGNs 


DIRTY 


MC 


Gordon+ 200 1 , Misselt+ 2001 


galaxies, nebulae 


TORUS 


MC 


Harries 2000, Harries+ 2004 


SE disks 


(no name) 


MC 


Heymann+ 20 1 2 


SE disks, AGNs 


SUNRISE 


MC 


Jonsson 2006, Jonsson+ 2010 


galaxies 


CRT 


MC 


Juvela+ 2003, Juvela 2005 


SE clouds 


(no name) 


MC 


Lucy 1999; 2005 


supernovae 


MCMax 


MC 


Min+ 2009; 201 1 


SE disks 


STSH 


MC 


Murakawa+ 2008 


SE disks 


MCTRANSF 


MC 


Niccolini+ 2003; 2006 


SE disks 


mcsim mpi 


MC 


Ohnaka+ 2006 


carbon stars 


MCFOST 


MC 


Pinte+ 2006 


SE disks 


HYPERION 


MC 


Robitaille 2011 


SE clouds 


PHAETHON 


MC 


StamateIIos+ 2004; 2005 


SE cores 


STEINRAY 


ED 


Steinacker+ 2003 


SE disks 




RayT 


Steinacker+ 2006 


SE cores 


(no name) 


ED 


StenhoIm+ 1991 


SE disks, AGNs 


HO-CHUNK 


MC 


Whitney+ 2002 


SE disks 


MC3D 


MC 


WoIf+ 2000, Wolf 2003b 


SE disks, SE cores, AGNs 


(no name) 


MC 


Wood+ 1999, Bjorkman+ 2001 


SE disks, galaxies 



^' finite differencing (FD), Monte Carlo (MC) or ray-tracing (RayT) 




Fig. 3. — Illustration of the various types of rays used 
in 3D dust RT calculations performed within an adap- 
tively refined density grid. Only a small fraction of the 
typical number of rays are shown for (a) an outer radia- 
tion source, (b) an inner radiation source, (c) a scatter- 
ing event, (d) a region with an optical depth near 1, (e) 
a coarse regular outer grid, and (f) rays to the observer 



pends on the specific needs of the application and/or 
the personal preferences of the user or developer. Nev- 
ertheless, the different approaches have some general 
strengths and weaknesses. These need to be consid- 
ered as rough guidelines only, as there are significant 
differences among codes based on the same approach. 
For example, simple MC codes based on the naive 
techniques explained in §5.1 are many orders of mag- 
nitude less efficient than modern MC codes that use 
weighting schemes. 

Generally speaking, setting up a reasonably effi- 
cient 3D MC RT code fits within the lifetime of a PhD 
project. Most MC RT codes are intrinsically 3D codes, 
and hence the shift from ID or 2D codes to a full 3D 
geometry is fairly straightforward. For RayT codes, on 
the other hand, the increase in complexity when mov- 
ing from ID and 2D to 3D is much steeper These 
differences explain the relative scarcity of general 3D 
RayT codes compared to MC codes in Table 1 ; a sim- 
ilar comparison of ID or 2D codes would result in a 
table with a much larger fraction of RayT codes. A big 
part of this complexity lies in the placements of the 
rays, which needs manual adjustment in RayT codes, 
but is done automatically in MC codes. The RayT pre- 
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Fig. 4. — Both simple and weighted Monte Carlo are 
illustrated graphically. This example includes 5 pho- 
tons resulting in one scattering photon reaching the 
observer and two absorption events for simple MC RT. 
This same set of 5 photons produces 5 scattered pho- 
tons towards the observer and 28 absorption events for 
weighted MC. The improvement in computational ef- 
ficiency of weighted Monte Carlo is clearly seen. 



calculation step needed for the manual placement of 
rays reveals an advantage for this solver as it identi- 
fies the critical locations in the model that dominate 
the appearance of the object. Other advantages of the 
MC method are that it needs less storage than RayT 
codes when scattering is included, and that it is widely 
used with many coders in the community improving 
its application to astrophysics with new algorithms, 
whereas RayT methods are mainly improved outside 
astrophysics. One strength of the RayT codes is the 
stronger and explicit error control. For example, the 
precalculation step can identify areas that are poten- 
tially under-resolved allowing for changes to the grid 
to provide adequte resolution for the full calculation. 
At this point, MHD codes tend to use RayT rather than 
MC techniques (e.g. Heinemann et al. 2006; Kuiper 
et al. 2012). 

7.2. Benchmark efforts 

Probably the most objective way to compare the 
strengths and weaknesses of the different codes is 
using benchmark problems. In the past few years, 
there have been benchmark efforts in many different 
computational astrophysics areas, including molecu- 
lar line transfer (van Zadelhoff et al. 2002), halo and 
void finder algorithms (Colberg et al. 2008; Knebe 
et al. 2011), astrophysical hydrodynamics (Agertz 
et al. 2007), cosmological hydrodynamical simula- 
tions (Frenk et al. 1999; O'Shea et al. 2005) and cos- 
mological radiation hydrodynamics (Iliev et al. 2006, 
2009). At the moment, no code validation or bench- 
mark project exists for 3D dust RT. The most advanced 
dust RT code validation projects are 2D benchmarking 
efforts. 

Pascucci et al. (2004) presented a benchmark test 
for two-dimensional equilibrium RT problems. Their 
system consists of a single star surrounded by a 
flared axisymmetric accretion disc. The optical depth 
through the disc varies from ry = 0.1 to Ty = 100. 
The accretion disc contains strong density gradients, 
which makes it an ambitious benchmark problem (un- 
fortunately, anisotropic scattering is not taken into ac- 
count). The authors compare the temperature maps 
and spectral energy distributions for five 2D dust RT 
codes (two grid-based codes and three MC codes). 
Differences between the various codes in the temper- 
ature maps are smaller than 1 % for the most optically 
thin model, but reach up to some 15% in the most 
optically thick system. For the emerging spectra, the 
differences range from a few percent for the optically 
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thin models to more than 20% for the most optically 
thick models. 

A first extension of the Pascucci et al. (2004) bench- 
mark test was presented by Pascucci et al. (2003)^, 
with two of the five codes from the original benchmark 
participating (one MC code and one RayT code, both 
intrinsically 3D codes). They considered a similar disc 
as in the Pascucci et al. (2004) benchmark, with an 
azimuthal ring added as a simplified model for a spi- 
ral density distortion. The main advancements were 
the addition anisotropic scattering and that images and 
visibilities, in addition to SEDs and temperature maps, 
were compared. The difference in flux in the entire 
image was smaller than 20%, but the visibility curves 
showed substantially larger discrepancies in their over- 
all shape. 

Pinte et al. (2009) go another step further, in what 
is the most advanced dust RT benchmark to date. They 
start from a similar circumstellar disc model, but con- 
sider optical depths up to Ty - 10^, use anisotropic 
scattering, and compare images and polarization maps. 
The four different codes that participate in the bench- 
mark are all MC codes. The agreement in the tempera- 
ture distributions is very good, with differences almost 
always smaller than 10%. Differences in the SEDs re- 
main smaller than 15% for models with Ty - 1000 
and agree within 20% over almost the entke wave- 
length range for the most optical thick cases. Pixel-to- 
pixel differences in high-resolution scattered light im- 
ages remain limited to 10% and the polarization maps 
do not differ by more than 5° in regions where the po- 
larization can be effectively measured by observations 
(Figure 7). 

8. THE FUTURE OF THE FIELD 

Present status 

3D dust RT is a rich and diverse field, with applica- 
tions across a broad range of astrophysical topics from 
dust near stars to entire galaxies. Correctly modeling 
the effects of dust on the transfer of radiation is criti- 
cal to studying many astrophysical objects, including 
the dust itself. Recent years have seen an impressive 
improvement in observational capabilities across the 
electromagnetic spectrum and this has shown that the 
dust distribution in many regions is strongly 3D. This 
requires methods to compute the dust radiative trans- 
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fer that can handle 3D structures and return solutions 
in a reasonable amount of time. The most common 3D 
dust RT solver is based on MC techniques, with RayT 
features in its modern accelerated form. A few appli- 
cations have used pure RayT solvers. Both methods 
face the challenges of grid discretization, determina- 
tion of uncertainties in solutions, and accurate com- 
parison between observations and the model calcula- 
tions. Almost 30 codes are currently able to deal with 
the full 3D dust RT with code variations arising from 
the prime field of application. There is no 3D dust RT 
benchmark, currently code comparisons are done us- 
ing 2D benchmarks. 

General trends 

Several trends indicate that the future of 3D dust 
RT is bright. The number of people actively involved 
in 3D dust RT is growing and the number of new 
published codes has increased significantly in recent 
years. A 3D approach to modeling complex distribu- 
tions is becoming common in many fields featuring 3D 
dust distributions. The continuing increase in available 
computing power and storage will support this trend, 
allowing a full transition from 2D to 3D dust RT for all 
objects showing 3D signatures. A prominent example 
of this trend is circumstellar disks with (proto-)planets, 
where the (M)HD simulations have been 3D for years, 
dust RT modeling often was 2D, and observations are 
now reaching the resolution to identify the 3D signa- 
tures of disk deformation due to a planet. In addition, 
modern online tools are expected to support the access 
to the codes by users through sophisticated interfaces. 

Future benchmarks 

For progress in 3D dust RT to continue, 3D dust RT 
benchmarks need to be established. Given the com- 
plexity of the codes, the increasing number of accel- 
eration algorithms, and large number of specific appli- 
cations, it is critical to provide a quantitative compari- 
son between codes. Experience with existing dust RT 
benchmarks and similar efforts in other areas indicate 
a suite of 3D dust RT benchmarks is needed. Ideally, 
each benchmark would focus on a particular part of the 
RT solution (e.g., scattering, polarization, equilibrium 
dust emission, or non-equilibrium dust emission) in a 
3D geometry. This would provide a clear test of differ- 
ent aspects of 3D dust RT and support the participation 
of all codes in at least part of the suite. 
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Data modeling future 

Given the impressive flow of new data now and 
coming years from ground- and space-based observa- 
tories, it is clear that the demand to accurately model 
3D dusty structures will rise strongly. Interfaces that 
can simulate observations with different telescope 
properties will become necessary to perform model- 
ing. We expect a rise of 3D dust RT modeling efforts 
that rely on automated fitting processes rather than 
by-hand explorations of the model parameter space. 
Since the number of sources of multi-wavelength data 
will rise, collaborations crossing borders between ob- 
servers and modelers will become more frequent. The 
utimate goal of 3D dust calculations is to modeling 
multi-wavelength images and derive quantitative and 
statistically sound information about 3D structures, 
embedded sources, and the dust itself. 

Future connections to non-dust RT codes 

Another future direction is the coupling of 3D dust 
RT codes with codes describing other physical effects 
in astrophysical objects. This trend is already happen- 
ing with 2D dust RT codes and the extension to 3D 
dust RT codes is clearly the next step. A variant of this 
type of connection is already happening where 3D dust 
RT is used to calculate the radiation field in a dust dis- 
tribution generated with a (M)HD code . Furthermore, 
(M)HD codes that make use of simple dust RT could 
be tested or the simple algorithms improved by com- 
parison to full 3D dust RT solutions. Chemical net- 
work calculations could be based on a more realistic 
estimate of the incoming radiation calculated from 3D 
dust codes. Finally, a combined calculation for 3D line 
and dust RT would enable line and continuum data to 
be simultaneously investigated using the same under- 
lying physical model. 

Future algorithms 

In the past, conferences and keyword related publi- 
cation searches have often been used to improve the 
unfortunately rare communication of new numerical 
algorithms from applied mathematics to astrophysics. 
The basic issue is the sheer flow of new findings and 
the different language of the two communities. Re- 
cent MC improvements have been mainly developed 
by coders working in the field, and additional efforts 
should be made to enable community-crossing ex- 
change on algorithms and error control. As a result 
of communications between coders preferring differ- 



ent solvers, we expect hybrid solvers making use of the 
advantages of the various approaches to appear more 
frequently. Given the increase in complexity in the 
modeled objects, we expect future activities to estab- 
lish grid generation algorithms that are optimized for 
3D dust RT; besides the octree or AMR-style grids that 
are now routinely implemented in 3D dust RT codes, 
unstructured grids as used in line RT (Paardekooper, 
Kruip & Icke 2010) and (M)HD codes (Springel 2010) 
are an interesting alternative. The inclusion of time 
dependence in the 3D RT problem, which could be 
important in the context of star formation or episodic 
accretion, will also need to be tackled with new algo- 
rithms (e.g. Harries 2011). The increasing availability 
of massively parallel machines will support algorithms 
that are optimized to run on many processors. 

Input physics improvements 

The improvement of the solvers is not restricted to 
developing algorithms that provide accelerated solu- 
tions. The interaction of radiation with cosmic dust 
is still not fully understood, and the variation of the 
dust properties with environment is an area of active 
research. The various continuum radiation sources like 
stars, PDRs, AGN accretion discs, and the ISRF are 
areas of vigorous investigation. For example, efforts 
based on existing and upcoming large scale surveys 
are being used to update the 3D structure of the stars in 
the Milky Way. Consequently, we expect to achieve a 
better understanding of the observed radiation from fu- 
ture research on the optical properties of dust, and im- 
proved data on the stellar and non-stellar sources that 
enter the 3D dust RT equation. 

Challenges 

A major challenge in 3D dust RT that emerges from 
this review is how to account for and mitigate sys- 
tematic uncertainties in the dust RT solution. They 
arise from under-resolving grids, from not propagat- 
ing rays/photons to important cells, and/or from uncer- 
tainties in the underlying dust grain models. As under- 
resolving of the dust and radiation field grid is often 
a result of constraints on computer memory and speed 
and improvements in algorithms to implicitly handle 
optimal grids are clearly needed. The pre-processing 
steps necessary for the RayT solver address some of 
these issues, but need further automating. The issue 
of not propagating enough rays/photons into particu- 
lar cells has been solved for both RayT (placement of 
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rays) and MC (biased emission), but in both cases cur- 
rently requires hand-tuning. An algorithm to automat- 
ically add additional rays/photons similar to that used 
for adaptive mesh refinement would clearly be useful. 
Finally, uncertainties in the assumed dust grain model 
provide a systematic uncertainty in the dust RT mod- 
eling that is difficult to quantify. Different dust grains 
models can be used to provide an estimate of this un- 
certainty, but the best way to reduce this uncertainty 
is to support the improvement of dust grain models 
through the use of improved laboratory and observa- 
tional data. 
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Fig. 5. — Left: Spitzer image of the molecular cloud 
L183 at 3.6 ;um revealing "coreshine" which is scat- 
tered light from the densest part of the cloud. Right: 
Spatial modeling based on basis functions which have 
Gaussian density structure in all three coordinates 
(Steinacker et al. 2010). The ellipses give the FWHM 
of the various Gaussians in the plane of sky. 




Fig. 6. — Examples illustrating the current state of the 
art in fitting dust RT models to images of galaxies. The 
left panels show the observed images, the correspond- 
ing panels on the right are the fits to these images. 
From top to bottom: (a) a clumpy 3D spiral galaxy 
model fit to an HST B-band image of the prototypi- 
cal edge-on spiral galaxy NGC891 by (Schechtman- 
Rook, Bershady & Wood 2012); (b) a 2D disc galaxy 
model fit to an SDSS g band image of NGC4565 by 
De Looze et al. (2012a) using a fully automatic fitting 
based on genetic algorithms (De Geyter et al. 2012); 
(c) a detailed 2D model for the Sombrero Galaxy, fit 
to a V-band image (De Looze et al. 2012b); (d) a 
3D clumpy disc model fit to an R-band image of 
the edge-on low surface brightness galaxy UGC7321 
(Matthews & Wood 2001). 
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Fig. 7.— An example of the Pinte et al. (2009) 2D dust 
RT benchmark, for a flared circumstellar disc with a 
V-band optical depth in the midplane of t = 10^. The 
image is a scattered light image, the three panels on 
the right show brightness profiles along the cuts plot- 
ted in the left panel, and the differences among the four 
different models used in the benchmark. The different 
line styles and colors indicate if the model was Pin- 
ball (Watson & Henney 2001), MCFROST (Pinte et al. 
2006), MCMax (Min et al. 2009), or TORUS (Harries 
et al. 2004). 
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